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Abstract 

We discuss heavy quarkonium production through parton fragmentation, including a 
review of arguments for the factorization of high-p-p particles into fragmentation functions 
for hadronic initial states. We investigate the further factorization of fragmentation func- 
tions in the NRQCD formalism, and argue that this requires a modification of NRQCD 
octet production matrix elements to include nonabelian phases, which makes them gauge 
invariant. We describe the calculation of uncanceled infrared divergences in fragmentation 
functions that must be factorized at NNLO, and verify that they are absorbed into the new, 
gauge invariant matrix elements. 
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1 Introduction 

The production of bound states of heavy quark pairs played an historic role in the development 
of the Standard Model p. This subject has retained a continuing fascination, in part because it 
offers unique perspectives into the formation of QCD bound states. The first step in quarkonium 
production, the inclusive creation of a pair of charm or bottom quarks, is an essentially per- 
turbative process. In particular, at high transverse momentum in hard-scattering processes, the 
dominant mechanism for the production of heavy quarkonium is the perturbative fragmentation 
of lighter partons, especially the gluon [2|. 
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A basic result of perturbative QCD for the production of hadrons H at high transverse 
momentum from the scattering of initial-state particles A and B is the factorization of universal 
fragmentation functions, [3 IH E] 

duA+B^H+xiPr) = da A +B->i+x(pT/z,fi) <g> D H/i (z,m c , //) + 0(m 2 H /p 2 T ) . (1) 

i 

Here, Cg> represents a convolution in the momentum fraction z. The cross section d<7A+B^i+x 
includes all information on the incoming state, including convolutions with parton distributions 
when A and B are hadrons, at factorization scale [i. We exhibit the m c -dependence of the 
fragmentation function Dn/i in anticipation of our interest in H as a charm- ant icharm bound 
state, for which the boundary condition for evolution is naturally taken at fi = 0(m c ). 

The actual transformation of a heavy quark pair into a heavy quarkonium H with momen- 
tum fraction z from parton i requires the introduction of fragmentation functions Dh/i{z, [i). 
This reasoning applies in principle to both charm and bottom quarks. For definiteness, we will 
generally refer to the heavy quark mass by m c , and to the produced hadrons as H = J/ip etc., 
of mass nin ~ 2m c . 

At moderate z only the evolution, that is, the /i-dependence, of the fragmentation func- 
tion is computable perturbatively, with the remaining information encoded in some initial func- 
tion D H / g (z, /i ), where we may take fi ~ ran- The effective field theory nonrelativistic QCD 
(NRQCD), however, has been invoked to simplify the nonperturbative content of the fragmenta- 
tion functions D H /i(z, Hq), in terms of a few (or anyway finite number of) nonperturbative matrix 
elements. 

The logic behind the application of NRQCD to fragmentation functions is readily summarized. 
One applies the NRQCD expansion in the relative velocity of the produced quark pair, assuming 
that a bound state will form only if this relative velocity is small to begin with. One then argues 
that the formation of the bound state is not affected by the exchange of soft gluons with other 
hard partons, only by exchanges between the quark and antiquark, and with the vacuum [7j. 
NRQCD then specifies a complementary factorization theorem, often written as 

daA+B^H+xM = ^ do A +B^cc[n]+x{PT) (0%) , (2) 
/( 

where the are universal NRQCD production operators, organized according to powers of the 
relative velocity of the cc state [n], and their rotational and color quantum numbers. We will 
encounter explicit forms below. 

Next, accepting that both (Q) and (j2J) hold at high-p^, the fragmentation function and 
NRQCD matrix elements are related by jHj 

D H/i (z,m c ,n) = ^2 di^ C c[ n }(z, //,m c ) (Cf ) , (3) 

n 

where dj_* C c[„] (z, n, m c ) describes the evolution of an off-shell parton into a quark pair in state [n], 
including logarithms of /i/m c . This formalism has been extensively applied to heavy quarkonium 
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production EH HU E21 Ed UH UH1 US]- At the same time, it has been observed that the 
applicability of NRQCD to production processes has not been fully demonstrated [To] IT7|. 

In this paper we revisit the formation of heavy quarkonia in fragmentation, with the aim of 
testing the relation Q. We will, however, first outline the proof of Eq. in leptonic annihilation 
and hadronic scattering, emphasizing that the arguments that apply to light quark bound states 
apply as well to heavy quarks, with corrections suppressed by powers of the mass of the quark 
divided by the transverse momentum. While our arguments for the factorization of fragmentation 
functions will not actually cover new ground, we are aware of no other unified presentation of 
the steps leading to (JTJ) for hadronic scattering in nonabelian gauge theories. 

Once we have established (0), we can test NRQCD factorization in its more specific form, 
Eq. (JHJ), which will simplify our analysis considerably. We study the factorization of the frag- 
mentation functions into perturbative coefficient functions times NRQCD matrix elements, with 
evolution included in the former. We shed new light on the necessary cancellation of infrared 
divergences in the perturbative calculation of coefficient functions. In particular, we will show 
that to carry out such a factorization, it is useful to redefine conventional NRQCD production 
matrix elements, with the addition of extra gauge links, or Wilson lines, a process that we termed 
"gauge completion" in Ref. |17| . This modification renders the matrix elements gauge invariant. 

Gauge completion is also necessary in order to absorb certain infrared divergences, beginning 
at next-to-next-to leading order (NNLO), that were not covered by the original arguments for 
NRQCD factorization. We should note that our NNLO infrared effects will appear only at order 
ctf in inclusive heavy quark pair production cross sections of the type calculated to order a 2 s in 
[T5] (see also the corresponding two loop decay cross sections in Indeed we will encounter 

two loop corrections with a quark pair and an additional hard gluon in the final state. The 
calculation is only possible, of course, because we restrict ourselves to the infrared sector. 

In any case, however, we are not yet able to prove Eq. Q to all orders in perturbation 
theory. The basic results of this paper were outlined in |17j . Here, in addition to the arguments 
on factorization, we will provide rather complete details on the methods used to identify the 
relevant infrared behavior, and on the necessary two-loop calculations. 

2 Factorization of Fragmentation Functions 
2.1 Long-distance dynamics in high-j9y production. 

The analysis of hard scattering cross sections begins with the identification of leading regions in 
the momentum space integrals of perturbative amplitudes and phase space 01201 E] Regions in 
this multidimensional space are conveniently classified in cut diagram notation, in which graphical 
contributions to the amplitude are represented to the left of the final state, and contributions to 
the complex conjugate amplitude to the right. In the complex conjugate graphs the roles of final 
and initial states are reversed. 

Because loop integrals are defined by contours in complex momentum space, it is only at 
momentum configurations where some subset of loop momenta are pinched that the contours 
are forced to or near mass-shell poles that correspond to long-distance behavior. These "pinch 
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surfaces" , or subspaces, in turn can be classified according to their reduced diagrams, found by 
contracting off-shell lines to points. 

The basic result is this. Any reduced diagram corresponding to a pinch surface can be 
interpreted as a physical process, in which each vertex can be assigned positions x% in space-time 
in such a way that if X{ and Xj are connected by line h carrying nonzero momentum ph, then 

Ax» = ( Xl - Xj f = ( Xl - xj) P£ , (4) 

where /3£ = p^/Ph is the four-velocity of particle h. At a pinch surface, this relation can be 
imposed for every line and vertex of the reduced diagram. Consistency then requires that the 
sum of Axjj's around any loop vanishes. This is enough to ensure that the reduced diagram does 
correspond to a physical picture, in which on-shell lines describe free, classical motion between 
the vertices. To this physical picture for finite-energy lines, lines with vanishing momenta may 
be attached in an arbitrary manner fZ0\ I2T] . 

The proof of the relationship between pinches in loop momentum space and physical pictures 
described by (jlj) is not difficult [201 1221 EHj , but we shall not review it here. Its consequences, 
however, are important and easily drawn. 

Let us apply the above considerations to the production of a hadron H, with momentum P, 
in the scattering of particles A and B, 

A(p A ) + B(p B )^H(P)+X. (5) 

We will assume that Pt is a large scale, of the order of the center-of-mass energy, and far above 
the strong coupling scale Aq C d- The relevant reduced diagrams for this process are illustrated 
in cut diagram notation by Fig. ^when A and B are a leptonic pair, and by Fig. |2] when A and 
B are strongly interacting (partons or physical hadrons). 




Figure 1: Reduced diagrams for high-pr particle production in leptonic annihilation. 
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At the pinch surfaces, there is a single hard-scattering, labelled by a shaded circle, in the 
amplitude and its complex conjugate. For dilepton annihilation, Fig. ^ the hard scattering 
is the result of the decay of the (real or virtual) electroweak boson formed in the annihilation 
process. For the hadronic process, the hard scattering results from the collisions of a single parton 
from each of the colliding hadrons 1 . Finite-energy partons emerge from the hard scatterings and 
connect to subdiagrams of on-shell collinear lines, the jets, Jj. All finite-energy particles of the 
final states are in one of these jets. In addition zero-momentum lines may be exchanged between 
the lines of the jets, and interact arbitrarily in a cut "soft subdiagram" S, consisting entirely of 
such lines. In effect, there are no final-state interactions involving finite momentum transfer in 
the reduced diagram for these processes. Thus, the total momentum of each jet is determined at 
the hard scattering, and the distributions of jet energies are calculable in perturbation theory. 
The essential observation to obtain this result is that once the parent partons of the jets emerge 
from the hard scattering, they separate at the speed of light, and subsequently cannot interact 
locally in any physical picture defined as above. The observed hadron H appears in one of the 
jets 2 . 

We will not attempt a full review of power-counting analysis in the neighborhood an arbitrary 
reduced diagram of Figs. ^ and El It is worthwhile to recall, however, that we may characterize 
these regions of momentum space by introducing scaling variables, conventionally denoted A, 
which control the relative rates at which components of loop momenta vanish near the pinch 
surfaces. In the terminology of [20], a leading region is one for which a vanishing region of 
loop momentum space near a pinch surface produces leading-power behavior. Such behavior is 
naturally associated with logarithmic singularities at the corresponding pinch surface. 

For the pinch surfaces at hand, we assign to the jth jet a light-like vector in the jet direction, 

lr There are, in fact, physical pictures corresponding to collisions involving more than one parton from each 
hadron. These, however, are suppressed by powers after summing over gauge- invariant sets of diagrams |24| . 
2 Here, we treat H as massless, on the scale of Pt- 
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n^, n 2 = 0, and an opposite-moving vector n^, n 2 = 0, normalized such that n,- • n,- = 1. For 
each jet, the combination of fij and rij define a two-dimensional transverse space, which we will 
denote as _L^. The fundamental leading regions are characterized by a familiar scaling behavior 
for the loop momenta internal to jet j, 

loop I in jet j : (nj ■ I ~ Ej, fij ■ I ~ XEj, l ± (j) ~ X l ^ 2 Ej) , (6) 

where is the energy characteristic of jet j, which we take to be of the order of the overall 
center-of-mass energy, denoted Q. In a similar notation, the scaling for soft loop momenta, 
internal to the soft subdiagram 5* or flowing between S and any of the jets, is 

soft loop : W~XQ. (7) 

In principle, the two scales A in Eq. (jUJ) and (J7j) need not be the same. A complete discussion 
includes power counting for subdivergences, as some lines approach the mass shell faster than 
others [201 HQ- 

2.2 Jet-soft factorization 

Arguments for the factorization of soft gluons from jets were given in some detail in Ref. [5] 
for e + e~ annihilation. We review these arguments here, and discuss their extension to hadronic 
scattering. Jet-soft factorization is made possible by the singularity structure of loop integrals 
near pinch surfaces jlj Many of these results have been rederived in the language of soft- 
collinear effective theory (211I2ZI- in graphical terms, the factorization is most clearly illustrated 
for leptonic annihilation, as in Fig. ^ 

2.2.1 The soft approximation in leptonic annihilation 

In Fig. 1 let us consider a loop momentum fcj, flowing from the soft subdiagram into jet Jj, 
through the hard scattering to another jet, and finally back to the soft subdiagram. We will 
examine poles in the integral near = due to the denominators of Jj. For definiteness, we 
assume the fcj loop is in the amplitude. 

The soft momentum appears in propagator denominators of a set of lines in the jet function. 
Let Tj be the momentum of any such line at ki = 0. When momenta are scaled as in Eqs. (jH)) 
and (|7j), any denominator in the jet function is of the general form 

(rj ± ki) 2 + ie = r| ± [ 2rij ■ rj fij ■ ki + Ifij ■ Tj rij ■ ki — 2rj± ■ ki± ] + k 2 + ie 
= 0(A) ± [ 0(A) + 0(X 2 ) + 0(X 3/2 ) ] + 0(X 2 ) + ie 
= r 2 ± 2nj ■ rj fij ■ ki + £>(A 3/2 ) + ie , (8) 

where the second equality gives the scaling behavior of each of the terms in order, according to 
(jHJ) and ((7j). In the third equality we exhibit the leading behavior, which depends only on the 
single component fij ■ ki of the soft momentum ki. This approximation will hold so long as the 
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contour of momentum component rij ■ k{ is not pinched in such a way as to violate the scaling of 
Eq. 0. 

Within jet Jj, the soft loop momentum ki can always be re-routed by shifts of the jet's loop 
momenta. In this way we may choose ki to flow on each jet line in a sense opposite to the 
direction of the jet momentum Pj oc fij (minus sign in Eq. (jHJ)). With this choice all singularities 
in the variable fij ■ ki are in the same (here upper) half-plane. Thus, although the poles in fij ■ ki 
due to jet lines are generally quite close to the origin, they do not pinch the integration contour. 
As a result, the presence of poles in the jet subdiagrams does not take us outside the scaling 
region of Eq. ((7j) . 

We now consider other possible sources of singularities in the variable fij ■ ki. Every virtual 
soft loop momentum that attaches at least one jet to the soft subdiagram may be routed so that 
it flows into no more than two jets. This is because once it reaches a second jet, it can be routed 
back to the first through the hard scattering subdiagram (where it is neglected in the off-shell 
propagators and vertices). The two jets determine distinct contour deformations for the soft 
loop momentum. These deformations are guaranteed to be compatible, however, because we can 
always specify them in a frame where the two jets in question, say j and j', are back-to-back. In 
this frame, we may identify fij ■ ki = kT and fiji ■ ki = kf ', for example. The consistency of the two 
contour deformations is then clear. Momentum ki must also flow through the soft subdiagram. 
Here, singularities are generically a distance 0(X) from the origin, except at lower-dimensional 
spaces. If such a subspace corresponds to a leading region, it may be treated separately, by 
the same arguments Finally, we note that if ki is the momentum of an on-shell gluon (or 
decays into a set of on-shell gluons), Eq. (jSJ) always holds, unless ki is itself collinear to the jet 
momentum. In this case, the line carrying ki should be treated as part of the jet. 

In summary, we have learned that the leading power behavior of the cross section from any 
leading region R may be found by keeping only the fij ■ ki component of soft momenta ki within 
the jet subdiagram, and setting the remaining components to zero. Similarly, in region R the 
soft gluons couple to the jet subdiagram only through the polarization component proportional 
to the jet direction, because at the pinch surface all other components of the jet tensor vanish as 
a power of A. 

Now we consider the subdiagram consisting of all lines in jet Jj, connected to the hard 
subdiagram by parton j in leading region R, not including the propagators of its external soft 
lines 3 . To be specific, we assume there are m soft gluons connected to the jet in the amplitude, 
and n — m in the complex conjugate amplitude. We introduce a jet function Jj R ' c ^ for leading 
region R, where C labels the particular cut of the jet subdiagram. With a given cut C, of course, 
the assignment of soft lines to the amplitude and its complex conjugate is specified. 

Our considerations lead us to a "soft approximation" for the function Jj R ' C ^ Within 
leading region R we may make a replacement that isolates the leading soft gluon momentum and 

3 In a covariant gauge, parton j is accompanied by a set of collinear vector lines with scalar polarizations at 
the coupling of jet J^' C ' to the hard scattering [5]. These gluons may be factored from the hard scattering and 
are included in the jet. 
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polarization components. In these terms, the soft approximation may be defined by 

J^ C) ZZl.. an (*i • • • *», ^i) = jf' C) ^a,..a n (*i • • • P) ^ • • ■ n« • • • n"" , (9) 

where we define, for any momentum ki, 

k~i = ki - fij rij (10) 

as a vector with only the "opposite moving" component of momentum. (Of course, the definition 
of ki varies from jet to jet.) The indices a* are the color (vector) indices of the external soft 
gluons of momentum ki, while M, L are the color indices of parent parton j, in the appropriate 
color representation. Corrections to the soft approximation are suppressed by powers of the 
scaling variable A, and hence by the overall hard scale Q. 

From the soft approximation Eq. (jUJ), the coupling of soft gluons to jet Jj is identical to the 
coupling of a set of unphysical gluons to the jet, whose polarizations are proportional to their 
momenta. Once we have made the soft approximation, it becomes straightforward to apply the 
nonabelian Ward identities of QCD to the connections of soft gluons to the jet [5]. This has a 
simple classical analogy. In the rest frames of particles within jet Jj, the classical fields due to 
particles in other jets, all of which are separating with relative velocities /3 re ; ~ c, reduce to pure 
gauge fields, up to corrections of order (/5 rei — c) [2H] • 



2.2.2 Factorization and the residual jet factor 



Once we have used the soft approximation and the Ward identities, the entire effect of the soft 
gluons external to the jets is to produce, order by order, a product of eikonal factors, 



1 {R,C)u 1 ...v n it, UP 
J j ML, ai ...a n I M • • • fWi> r 3 



n vi ... n Vn n^ 1 ... 



j( R > C )(p\ pO') Mm+l-Mn t I _U 
J j \ r 3) rj MK,a m+1 ...a n \ 



- K ) ( kl 



kr. 



Here, m is the number of soft gluons that couple to the jet subdiagram in the amplidue, and 
n — m the number in the complex conjugate amplitude. The eikonal factors E and reproduce 
all momentum and color dependence, but are insensitive to the internal dynamics of the jet, and 
depend only on the 4- velocity fij, in the jet direction, and the color representation of parton j. 
Specifically, they are given by 



p(j) HI. ../J, n 



k\ . . . k^ 



E 



perms 



ghfT, 



i=l 



-n~ 



(k 1 + --- + k i )+ ie 



(12) 



KI 



where P implies ordering of the color matrices T'i"'' according to the permutation of soft gluon 
connections. (As above, soft momenta flow into the jet.) The function Jj R ' C \Pj) in (jlljl repre- 
sents what we will refer to as the residual jet factor in region R with cut C. It is given by the 
normalized color trace of the jet function with no external soft gluons, 

1 



J 



(R,C) 



(Pi, 



d(j) 



(R,C) 



LL 



(Pi) 



(13) 
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where d(j) is the dimension of the color representation of parton j. 

To apply the Ward identities that lead to Eq. (|TT|) we need only integrate over the opposite- 
moving components, nj ■ I of the jet loop momenta I. This is because the Ward identities only 
require shifts in loop momentum equal to the momenta that flow into the jet from the external 
lines. 

In general, the residual jet function includes contributions from soft gluons for which the 
soft approximation fails, but which remain internal to the jet. It is not necessary that the soft 
approximation apply to every soft gluon. Rather, for this analysis to hold it is only necessary 
that in every leading region we can find a set of soft lines for which it holds, and for which we may 
apply Eq. (JTTJ). Equation (JTTj) is a general result for final-state jets in arbitrary leading regions. 
We will see below how function J- ' can be identified as a contribution to a fragmentation 
function. 

At each (here nth) order, the factorized hi dependence in the eikonal factor E in (|12|) is 
identical to the corresponding dependence in the expansion of the ordered exponential 



<Dg(0,oo) = Pexp 



ig / dXfij ■ A (j) {n\) 



(14) 



where now P denotes path ordering, and where A^> is the gauge field in the matrix color repre- 
sentation of the parent parton of the jet (quark, antiquark or gluon) 4 . At leading power in A, 
and hence in the large momentum scale Q, soft gluons couple to jets only through the operators 
Uj ■ A, restricted to the light cone along the jet directions. We will see below how other operators 
arise at nonleading powers. 

All of the reasoning above may be applied to the particular jet (Ji in Fig. 1) from which 
the observed hadron H arises. The entire leading-power dependence on the masses and relative 
momentum of the quarks, as well as on the momentum fraction (z) of the pair is in the functions 

( Ft 

J- ' (P) at each leading region. The influence of soft gluon emission on z can be neglected, 
precisely because of the soft approximation Q. Thus, in each leading region, the jet dynamics 
that produces an observed particle decouples from soft gluons that could link it to the other jets 
in the final state. In this way, fragmentation is seen to be universal, depending only on the parent 
parton, the produced hadron, its momentum fraction z, and eventually a factorization scale. 



2.2.3 Hadronic scattering 

The arguments for jet-soft factorization in hadronic scattering are similar to those for leptonic 
annihilation 5 , but special care must be taken because of the "initial-state" jet subdiagrams J a 
and Jb consisting of lines collinear to particles A and B in Fig. EJ As we shall see, however, the 
factorization of fragmentation within a final-state jet holds in this case as well, and is actually 
somewhat more general than collinear factorization in terms of parton distributions jHj- 

The essential difference between leptonic annihilation and hadronic scattering may be seen 
by comparing Figs. ^ and |21 In the former, although the poles from jet subdiagram Jj in the 

4 This is the factorization effected in soft-collinear effective theory by a redefinition for collinear fields. 
5 They are also essentially identical to those for lepton- hadron scattering. 
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soft momentum components fij ■ ki are closer to the origin than 0(A) in general, they are all in 
the same half-plane. As a result, these momentum components may be deformed away from the 
poles into a region where the soft approximation holds. 

For hadronic scattering, precisely the same reasoning applies for soft momenta that flow only 
between final-state jets, and/or through the hard scattering. It also applies for soft loops that 
connect to an initial-state jet only via lines whose large momenta flow directly from the initial 
state into the hard interaction. In these connections, to what are sometimes called "active" jet 
lines, all poles are again in the same half-plane, and the same reasoning allows us to deform 
contours as above to justify the soft approximation. 

A difference arises, however, when soft lines connect to the initial state jets by "spectator" jet 
lines, whose momenta flow into the final state without passing through the hard scattering. In 
this case, to complete the soft loop through the hard interaction, the momentum must flow "back" 
to a vertex at which spectator lines and active lines connect, and then flow once again forward 
into the hard scattering. Suppose this occurs for initial-state jet A. Both spectator and active 
lines in A produce poles close to the origin for a soft component ua- ki, and these poles are in 
opposite half-planes. The resulting pinch forces us into a leading region where Ua • ki -C k i± ~ A, 
which is generally referred to as a "Glauber region" [311221 • in this region, the scaling (jjj) does 
not hold, and the soft approximation fails for this jet. Because the soft approximation fails, soft 
gluons "resolve" the internal structure of the jet, and the factorization arguments given above 
may not apply. 

When the soft loop flows between an initial state jet and a final state jet, however, only a single 
light-cone component is pinched, associated with the initial-state jet. The soft approximation 
may still be applied to the final-state jet, giving eikonal factors as in Eqs. (J 11)) and The 
eikonal factors associated with the outgoing jet then cancel in a single-particle inclusive cross 
section, in same way that soft divergences cancel in jet cross sections. 6 This decoupling and 
cancellation of soft gluons enables us to identify universal fragmentation functions, in terms of 
universal matrix elements, in hadronic scattering as well as leptonic annihilation, independent of 
the jet structure of the particular hard scattering. 

Finally, consider those soft loops that flow between the two initial-state jets and/or through 
the hard scattering. In general such loops encounter Glauber pinches in two light-cone compo- 
nents. For cross sections that are inclusive in soft gluon emission and in the fragmentation of 
the forward jet remnants, these pinches nevertheless cancel in the sum over final states. We then 
have collinear factorization into independent parton distributions for incoming hadrons A and 
B and a fragmentation function for hadron H It is worthwhile noting, however, that even 
when these criteria are not satisfied, and the overall cross section does not factorize into incom- 
ing parton distributions (as, for example, in diffractive scattering in hadron-hadron collisions 
[3T]). the final-state jets still factorize from the incoming jets and their soft exchanges, and the 
single-particle cross section at high is still governed by a universal fragmentation function. 

6 The argument for this cancellation in the case of hadronic scattering is given in the first part of Sec. V of Ref. 

POl- 
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2.3 Power corrections 



Once we have determined that the leading-power contributions factorize for the leading regions 
associated with Fig. |21 we naturally turn our attention to power corrections These may be 
classified by an expansion in nonleading contributions to the integrand near the pinch surfaces. 
It is therefore an expansion in terms of ratios such as f(k)/qj(k), where the numerator f(k) 
represents any of the terms involving soft momenta k that can be neglected at leading power. 
These are the terms in Eq. (|SJ) that scale as A 3 / 2 or higher, as well as non-leading terms from 
numerator momenta. The denominator represents the squared momentum of the jet line qj, after 
the soft momenta k flowing on jet line qj is replaced by k = nj(fij ■ k) ~ 0(A), Eq. (fTUj). Although 
q"j{k) is not large, the ratio f{k)/q? is small in the leading region. (As we have seen, this may 
require contour deformations.) 

Keeping only nj ■ k terms in the denominators, the numerator terms f(k) are polynomials in 
the rij ■ k and components of soft gluon momenta. In position space, these vertices, connected 
to jet lines, correspond to operators that are local with respect to the rij and directions, 
but are relatively on the fij light-cone. The gauge invariance of the theory requires that these 
vertices, representing the interactions of soft gluons with the jet functions, combine to form 
gauge covariant operators. We may think of these vertices as supplementing the leading-power 
fij ■ A vertices of the "soft approximation", identified above with the Wilson lines of Eq. (|14p. As 
above, the application of Ward identities, or equivalently a redefinition of collinear fields as in 
soft-collinear effective theory, organizes all leading vertices into nonabelian phase operators, but 
now acting as color rotations on the nonleading vertices as well as on the "parent" parton lines 
of the jet functions. For the purposes of factorization at leading power in Q, however, we need 
not enumerate these nonleading operators or vertices. 

2.4 Fragmentation functions 

So far, we have identified the leading regions in cut diagrams that are associated with infrared 
dynamics in single-particle inclusive cross sections. We have seen that at each leading region 
the cross section breaks up into a factor associated with the production of a parton (j above), 
times a jet function Jj R ' c \P) that describes a contribution to the formation of hadron H. In 
this section, we will show that the jet functions identified above are in one-to-one correspondence 
with leading regions for the standard fragmentation functions, Dh/j- 

Fragmentation functions may be defined in terms of expectation values [33]. For example, 
consider the production of hadron H from a parent gluon at momentum fraction z, taken, for 
definiteness along the 3-direction. The relevant matrix element is then, in D = 4 — 2e dimensions, 

D H/9 (z,v) = 16{2 ~^l)nP+ rRcolOT / dx ~ e ~ iiP+/Z)X ~ (0]^ +A (0) [^ A) (0,°o)] f 

xa H (P + , ± )4(P + , ± )$i A) (0, oo) F + A (0+, x~, 0±) |0> , (15) 

where a} H is the creation operator for particle H at momentum P and F +x is the gluon field- 
strength. The operator is defined as in Eq. (|T4*j) . but in the direction = 5^-, opposite 
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to the direction of hadron H. Its fields are in the adjoint matrix representation of color. The 
product of operators on the light cone requires renormalization and the introduction of a scale 
[A, as described in [53] . 

The expectation value in Eq. (jl5j) may also be expressed as a sum over all states including 
hadron H, 



This form shows its close correspondence to a cross section. 

The leading regions of the expectation values (fTK|) and (fTT)|) are, in fact, very similar to those of 
leptonic annihilation cross sections discussed above. Every leading region includes in its reduced 
diagram a jet J{P) that provides the particle of momentum P in the final state, in addition to 
a jet in the opposite-moving direction n^, and possibly other jets and arbitrary soft radiation 
(subject to the effective phase space limitations imposed by renormalization at scale /i). At each 
such leading region R, the same arguments as for leptonic annihilation lead to the precise analog 
of Eq. (JUJ), with exactly the same residual jet functions Jj R ' C \p). Because the cross section is 
otherwise inclusive, the sum over final states results in the cancellation of all soft and collinear 
singularities except for those associated with J- ' (P). Therefore, we recognize a one-to-one 
matching of every leading region in the fragmentation function with a corresponding region in 
the total cross section. This is the case for both leptonic and hadronic initial states, because the 
residual jet functions are the same in each case. 

Strictly speaking, of course, the above discussion applies only to perturbation theory, which 
requires that we impose an infrared regulation, presumably dimensional regularization. Because 
our arguments extend to all orders in perturbation theory, however, we may in principle intro- 
duce an interpolating field with the quantum numbers of hadron H, sum to all orders in that 
channel, and isolate the S-matrix elements for H in the regulated theory. In this sense our ar- 
guments demonstrate factorization for bound state H in the regulated theory. We assume that 
the continuation back to physical QCD in four dimensions respects this result. This assumption 
is shared with essentially all demonstrations of infrared safety and factorization. 

3 NRQCD Factorization and Gauge Completion 

Having reviewed arguments for the factorization of fragmentation functions, Eq. up to cor- 
rections in powers of tuh/Pt, we are ready to rephrase the question of NRQCD factorization in 
terms of the fragmentation functions themselves, as in Eq. ©. We begin with a further examina- 
tion of the leading regions of the fragmentation functions, and we discuss evolution to the mass 
scale of the heavy quarkonium m#. We then analyze the refactorization, Eq. Q of the gluon 
fragmentation function in terms of NRQCD production operators, and propose a gauge-invariant 
extension of the conventional operators. 




N 



x | H(P + , ± ) N)(NH(P + , ± ) \& n A \0, oo) F+ A (0+, x~,0 ± ) |0> . (16) 
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3.1 Refactorization at the heavy quark mass 

Our first goal is to separate logarithms associated with evolution from dynamics at the scale 
of the heavy quark mass. This can be done by invoking the evolution equations for the gluon 
fragmentation functions in Eqs. (fTSJl and (JTfiJl . 



with a sum over partons z, and similarly when the gluon is replaced by a quark or antiquark. 7 
The solution to ()17|) enables us to relate fragmentation at the conventional scale Pt with the 
mass scale of the produced hadron, H, 



where C gi is a perturbative factor. 

We will want to study the expansion in relative velocity of the heavy quarks in the fragmen- 
tation function evaluated at a scale on the order of the heavy quarkonium masss. It is natural, 
of course, to carry out this expansion in the rest frame of the heavy quark pair. Since this is 
not the usual frame in which to discuss fragmentation or the evolution (|17|) associated with it, 
we will briefly discuss how evolution appears in this frame. Specifically, we need to show that 
evolution logarithms factorize from the decay of an off-shell gluon, with mass of order m#, as 
seen in the rest frame of hadron H. 

The transverse momentum of the observed heavy quarkonium in the fragmentation function 
(J 15)) is by definition zero. Thus, the transformation to its rest frame is a boost in the direction of 
its momentum as seen in the lab. For convenience we take this momentum in the "plus" direction, 
as in (|15|). In both the lab frame and the quarkonium rest frame, evolution then results from the 
strongly ordered transverse momenta of partonic radiation. 

To confirm Eq. (|18|) . we should verify that we can factorize soft gluons that connect partons 
with transverse momenta kj_ 3> vrtu from those of lower transverse momentum. The former will 
appear in the evolution functions C g i, the latter in the fragmentation function at the scale of m#. 
This separation of low- from high-fc^ gluons as seen in the H rest frame follows exactly the same 
pattern as the factorization of gluons from the jets in Sec. 2 above. 

Consider a parton d of transverse momentum k^± 3> m# and longitudinal momentum fcj" = 
ZdQ, as seen in the lab frame (or the center of mass frame of the overall collision), with Q the 
energy of the jet in that frame. In the same frame and notation, the heavy quarkonium H has 
transverse momentum k±^ = 0, and energy Eh = zQ ^> run- A boost to the rest frame, where 
the energy of H is m#, leaves the transverse momentum k^± unchanged, while transforming the 

7 The evolution kernels for heavy quarks may be chosen identical to those for massless quarks in the case of 
parton distibutions [H4| . A similar relation should hold here, although we will not attempt a formal proof. 




(17) 
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plus and minus components of k d according to 



k~ = ^ => ^ (19) 
d 2z d Q 2z d Q m H y ' 



Equivalently, the rapidity of parton d transforms according to 



By assumption, k d ± ^> rriH- Therefore, as long as z is not itself small, that is, assuming that H 
is one of the "leading" hadrons in the jet, the rapidity of parton d, which is large and positive in 
the center of mass frame, is large and negative in the rest frame of hadron H. In this frame, all 
strongly ordered (in transverse momentum) partons are moving in the direction opposite to the 
original jet direction. The soft approximation can now be applied to soft gluons connecting the 
heavy quark pair that forms the quarkonium to the strongly ordered gluons. The interactions of 
these soft gluons may then be approximated by an eikonal line in the direction n M , opposite to 
the jet's direction. The only difference from soft-jet factorization in a cross section is that now 
the soft gluons' transverse momenta are smaller than m#. The result is exactly a fragmentation 
function with upper limit mn on gluon transverse momentum in convolution with a perturbative 
function, as in Eq. (fTH|) . which is what we set out to show. 



3.2 Long and short distance dependence at the scale ra# 

To make contact with NRQCD applied to a fragmentation function, we explore further the sources 
of its long- and short- distance behavior. This can be done as in the discussion of cross sections 
and fragmentation functions above, in Sec. 2, although now we will carry out our analysis in the 
rest frame of the heavy quarkonium. We begin, as above, with the physical pictures associated 
with pinch surfaces. 

The relevant physical pictures for fragmentation into hadron H are shown in Fig. El Since 
we are working in infrared regularized perturbation theory, the heavy quarks appear in the final 
states. We recall our discussion above, however, in which we argued that in principle the reduction 
of the bound-state pole does not modify factorization. We will continue with this assumption. 

A related point is that at the bound state pole the relative momenta of the quark pairs on 
either side of the cut need not be the same. In principle, then, we should take the relative 
momentum of cc pair in the amplitude, q below, to be independent of the relative momentum, q' 
of the cc! pair to the right. This is the method employed in the explicit calculations of Refs. jHSJ 
EH! E3 EEl E3 HO] > for example. Powers of q and q', however, are employed to identify operators 
in NRQCD, terms linear in q corresponding to the lowest order of the covariant derivative. Since 
we are interested primarily in separating infrared poles from coefficient functions, we will not 
distinguish between q and q' below, and simply calculate the fixed-order eikonal cross section for 
a quark pair. 
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Figure 3: Physical pictures for heavy quarkonium production. 



As noted in the discussion of Sec. 2.4, the physical pictures for fragmentation are similar to 
those for the hadronic final state interactions in leptonic annihilation. The process begins with a 
short-distance subdiagram, represented by a shaded circle in Fig. El In this case, short-distance 
refers to virtualities at the order of m#. Lightlike jets, Jj, develop from energetic (E ~ tuh) 
"semihard" quarks, antiquarks or gluons, which emerge from the hard scattering. These jets may 
be connected by a diagram consisting entirely of soft quanta, S, to each other, and to the Wilson 
line that is part of the construction of the fragmentation function. In the case at hand, the heavy 
quark pair also emerges from the hard scattering, and soft quanta may connect to the jets and/or 
the Wilson line. To prove NRQCD factorization, Eq. Q, it will be necessary to show that all of 
this long-distance behavior either cancels or matches entirely to NRQCD matrix elements. 

Considered abstractly, the connection to NRQCD is made by "integrating out" degrees of 
freedom at the mass scale m# in the calculation of the fragmentation function. In practice, 
that is in perturbation theory, the NRQCD operators O n can be identified once we consistently 
separate long and short distance contributions. As the figure shows, a generic pinch surface in 
phase space involves not only a truly short distance part, but also a variety of semi-hard jets. 
The question we must ask is to what extent hadronization is affected by the presence of these 
jets. In the original discussion of NRQCD factorization given in Ref. j7], it was argued that in the 
inclusive sum over cuts in H production, all infrared divergences due to soft exchanges between 
the heavy quarks and the extra jets cancel in the inclusive sum, even while we fix the final state 
of the quark pair to be a gauge singlet. Notice that even in the absence of semi-hard gluons, soft 
gluons may be exchanged with the Wilson line that is part of the definition of the fragmentation 
function. Indeed, this Wilson line is what remains of all exchanges of soft gluons between the 
heavy quarks and partons at relative momenta greater than m#. 

In the absence of the soft gluon connections between the heavy quarks and semi-hard gluons, 
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the remaining physical pictures can be represented as in Fig. In this case, all collinear and soft 
divergences associated with the jets, whose final states are summed over inclusively, cancel, just 
as in leptonic annihilation. Soft singularities may, and in general do, remain in the transition of 
the heavy quarks from short distances to hadronization, but such soft divergences are said to be 
"topologically factorized" [121 > anc ^ are readily factorized from the hard scattering function by a 
standard expansion in relative velocity, as we now sketch. 




Figure 4: Topologically factored physical pictures for heavy quarkonium production. 



In a topologically factorized diagram like Fig. we can expand the short- distance function 
around vanishing relative momentum, or equivalently relative velocity v of the heavy quarks. 
Similarly, we may decompose each diagram according to the color state (singlet or octet) of the 
heavy quark pair, and may also expand in the momenta of any light quanta (gluons or quark 
pairs) that also emerge from the short distance subdiagram. This leads precisely to an expansion 
in terms of local operators, creating heavy quark pairs in states n, labelled cc[n] , where in general 
[n] also labels the term in the expansion in relative velocity and light parton quanta, so that the 
corresponding operator, ip^ (0) n' n x(0) , which always includes a quark pair, may also create light 
quanta. Each such operator will be accompanied by the sum of all hard subdiagrams, evaluated 
at zero relative velocity and at zero light parton momentum. We will refer to this sum as the 
hard scattering, or coefficient, function for operator ip^ (0) n' n x(0) . 

Combining the expansions from the amplitude and its complex conjugate, we derive Eq. (J3J), 
with operators that describe the creation of a heavy quark pair from the vacuum, summing over 
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all final states that include hadron H . The general form of these operators [7j is 



N 



= ^(0)^(0) (a^o H ) ^(0Kx(0) , 



(21) 



where the insertion of the creation operator a H , which is understood to act on out states to 
produce hadron H, and its conjugate enable us to sum over the complete set of out states 
between the creation and annihilation operators. The first form defines the sum over final states 
appropriate to quarkonium production, while the second form is a convenient shorthand. 

The following discussion is an attempt to analyze the basic assumption that enables us to 
expand in v in this manner. That is, we will begin with the general momentum region illustrated 
by Fig. El and explore the reduction to the simpler "topologically factorized" picture of Fig. HI 
by testing the cancellation of soft exchanges between the heavy quarks and semi-hard gluons, or 
equivalently, the Wilson line. 

3.3 Operators and gauge completion 

Our first observation, already described in is that matrix elements of the form (|21jl are not 
invariant under operator-valued gauge transformations. In general, the onium creation operators 
an and a* H , which act on out states, need not commute with gauge transformations carried out 
at the origin, even though they are themselves color singlets. As a result, it seems most natural 
to us to modify the operators (J2T]) to provide a form precisely analogous to the gauge-invariant 
definitions of fragmentation functions in Eq. (|15J1 above, 



in terms of ordered exponentials, defined as in Eq. (|14|). In the (complex conjugate) amplitudes, 
(anti)time-ordering is understood. We emphasize that such a redefinition is not required for 
self-consistency. If one can demonstrate NRQCD factorization in terms of operators in any 
specific gauge, a gauge-dependent definition of the operator matrix elements is admissible, as 
long as the gauge-dependence is not infrared sensitive. Indeed, this is the case for fragmentation 
functions, because of the cancellation of infrared divergences in final-state interactions at high 
Pt-, as observed above. In the absence of a similar demonstration of infrared finiteness for the 
refactorization Q of fragmentation functions in terms of NRQCD operators, however, it seems 
natural to entertain (J22|) as a plausible replacement. We now turn to the expansion in relative 
velocity, which will enable us to test our suggestion. 



(22) 
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4 Velocity expansion 



4.1 Requirements for NRQCD factorization 

To study the role of soft gluon emission in heavy quarkonium production, we will analyze infrared 
divergences in the production amplitude for two heavy quarks, of total momentum P and relative 
momentum, q: 

P P 
P 1 = — + q = p + q p 2 = — - q = p- q . (23) 

That is, we study the process g — > cc[n Q ] + X with cc[n ] = c(p + q)c(p — q). The lowest-order 
diagram for this fragmentation function is shown in Fig. |S| It consists of a single gluon splitting 
into the quark- ant i- quark pair. 



P/2+q 




Figure 5: Lowest-order fragmentation function for g — > cc. There are no interactions on the 
eikonal quark pair or the Wilson line that corresponds to an eikonal gluon of four-velocity I. 



In the following, we will study infrared divergences in soft gluon corrections to this process, 
when the quark-antiquark pair is created as a color octet, but is restricted to a singlet in the 
final state. Otherwise, we sum over all perturbative final states. 

To form a heavy quarkonium, of course, these quarks cannot be truly on-shell. Rather, they 
are off-shell by an energy of order q 2 /m c , characteristic of a Coulomb bound state [H]. These 
are nonperturbative effects, however, while coefficient functions are calculated in perturbation 
theory. The cancellation of divergences, and/or their matching to matrix elements in soft-gluon 
corrections to Fig. El is a necessary condition for NRQCD factorization. Any remaining diver- 
gences would be a violation of factorization. In our calculation below, we will find uncanceled 
divergences at NNLO for conventional operators, which, however, may be absorbed into gauge- 
completed NRQCD operators. 

4.2 Expansion in the eikonal approximation 

Because our calculation will be carried out with on-shell quarks, we can use the eikonal approx- 
imation for the coupling of soft gluons to the quarks in order to identify infrared divergences 
in the cross section. Equivalently, we may treat the quarks in heavy-quark effective theory to 
leading order in their mass. Yet another equivalent, and for us particularly convenient, approach 
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is to replace the quarks by path-ordered exponentials, similar to Eq. (|14j) above, but now with 
time-like velocities representing the quark and antiquark. 

The dimensions of the velocity in ordered exponential (fTlj) can be shifted by a change of 
variables in the parameter A. For this reason, we are free to identify the quark velocities directly 
with their momenta = (P/2 ± q). At fixed (and unequal) values of Pi and P2, all infrared 
divergences can be found by the eikonal approximation. The eikonal approximation and hence 
infrared divergences are completely independent of any spin projections that we may make on 
the state of the quark pair. As a result, soft gluon emission separates from all dynamical factors 
that involve the spins of the quarks, and enters as a multiplicative factor. We will come back to 
the limit Pi = P2 below. 

We can classify the eikonal infrared-sensitive factors by the color of the cc pair at creation 
(the origin of the ordered exponentials) and by their color in the final state. Gluon emission, of 
course, will mix these states. For an NRQCD-like factorization to hold, if we fix the color of the cc 
pair in the final state, infrared divergences either cancel or can be matched with matrix elements 
[7j. Finite remainders will be associated with coefficient functions, as in the NLO calculations of 

Refs. Esj-gn] 

In summary, we will study the infrared factor associated with the creation of a cc pair in an 
octet configuration, and its evolution into a singlet in the final state. This infrared factor may 
be written in the notation of Eq. ()14|) as 

J^(Pi,P 2 ) = £<0| \^(0)]jT d ] JK [$gt(o)j ^ ^(0) db \N) 



N 



x(N\^ A \0) bc *w (0) [T C ] MN <&«,'(()) ^ |0>, (24) 

where we have exhibited all color indices: those in adjoint representation by a, b ... , and those 
in the fundamental representation by I, J . . . , to indicate the trace structure, which imposes a 
color singlet configuration in the final state. 

The operator $® is the ordered exponential that represents the antiquark. It has the opposite 
sign on the coupling compared to the quark operator, and has color matrices ordered in the reverse 
sense to time ordering. In the notation of the standard definition, Eq. (|14j). we represent this 
matrix ordering by V, and define 



Kl 



■M, 



*g(0) 



Pexp 



d\P 2 - A iq \P 2 \) 



(25) 



Here = T a A v> a is the matrix- valued field in the quark fundamental representation. For 
classical fields, $n (0) is the hermitian conjugate of $1^(0). In Eq. (J24|) and below, overall 
time-ordering of the field operators is understood in the amplitude, and anti-time ordering in its 
complex conjugate. For explicit computations, we restrict the sum over final states iV in Eq. (}2"4"j) 
to soft gluon emission only. 

The graphical rules for the interactions of gluons with the ordered exponentials are exactly 
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the same as the eikonal approximation, and propagators and vertices are given by 



(j3-k + ie) 



±ig s T a (3» , 



(26) 



with the plus for antiquarks and the minus for quarks on the vertex and with (3^ the time-like 
quark four-velocity. The quark and antiquark eikonal propagators are represented as heavy lines 
on the left-hand side of Fig. |H1 In this notation, Eq. (jUj) describes a product of color traces in the 
fundamental representation. Our ability to use the same notation for velocities as for momenta 
is manifest since the combination of each eikonal vertex and propagator is scale- invariant. In 
the pair rest frame the relative velocity of the members of the pair is proportional to the ratio 
q 2 /m 2 c = Aq 2 /P 2 . 

In the spirit of NRQCD analysis, and because it leads to some simplification, we will study 
corrections to Fig. El to order q 2 , which is the first nontrivial order. At zeroth order in q 2 , the 
quark and antiquark never separate, and all infrared divergences cancel, since there are no color 
multipoles to which they can couple. We can see this in Eq. ()24jl . in which both the amplitude 
and complex conjugate amplitude reduce to unity in the limit P\ = P2 — ■> P/2. This is easily 
proved by considering the A(x)-field with the largest time in the amplitude. This field may come 



either from the quark exponential, $^(0), or the antiquark exponential $p 2 j (0). When Pi = P 2 , 
the only difference between these two terms is the relative minus sign between the quark and 
antiquark vertices. Every such pair of terms cancels pairwise. An identical argument applies to 
the complex conjugate amplitude, and there is therefore no overall q term, and q 2 can be reached 
only by expanding both the amplitude and its complex conjugate to order q independently. 

The expansion to order q is straightforward, and has a nice interpretation in terms of fields. 
We start with the expansion for the individual ordered exponentials, 
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(27) 



where the explicit minus sign on the left in the second expression anticipates that we will be 
expanding in the momentum of the anti-quark, P2 = p — q. In both of these expressions, the 
time-ordering is from right (earlier) to left (later), with an (opposite) identical ordering of color 
matrices for the (anti) quark exponential. We have inserted an explicit T in the antiquark 
expression, to remind ourselves that the operators and color matrices have the opposite ordering 
in this case. The operator F^ VA is the gluon field strength with tensor and color indices. Note 
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the overall factor of A', which reflects the increasing separation of the quark and antiquark paths 
with increasing distance from the origin when q is changed by a constant amount. 




+ 





+ 0(v2) 



Figure 6: Velocity expansion. The heavy solid lines represent heavy quark propagators in the 
eikonal approximation. 



We now apply Eq. (|27jl to the amplitudes in Eq. 
we find 



Expanding Pi and P 2 about p = P/2 
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NL 



Here again, time ordering is understood for all field operators. The lowest order of the expansions 
of the left- and right-hand sides of Eq. (|2*K|) are shown graphically in Fig. El On the right, the 
vertex associated with the field strength in the final form is represented by ®. 

We next apply reasoning similar to that which led to the cancellation of the ordered expo- 
nentials at q = 0. Again, consider the A-field with largest variable A, assuming that there is at 
least one such field with A > A', that is, at least one field at a larger time than the field strength 
p il q v F vpi . We recognize that whenever we find such a field, there is a cancellation between the 
cases when that field is associated with the quark and antiquark ordered exponentials. All fields 
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at times greater than that of the field strength cancel, and we have 
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In the last line we have used the cyclic nature of the color trace and the anti-path ordering of 
the antiquark exponential to show that the two terms above are equal. Next, we note that were 
it not for the generator T c , we could use the same reasoning as above to show that the A-field of 
lowest A' cancels between the quark and antiquark exponentials. The nonvanishing remainder, 
therefore, is a commutator, 



-igA d (X min p) [T d ,T c 



9 Ad (X min p) fdceXe 

-igA d (X min p) [T { d A) } ec T e , 
T e (~zgA^(X mmP )) ec , 



(30) 



where [Tj^] ec = —if dec is a generator in the adjoint representation. In effect, the gluon field is 
converted from the fundamental representation to the adjoint. 

At any order in g, this procedure may be repeated until all A-fields from the remaining 
ordered exponentials have been converted from fundamental to adjoint representation in Eq. 
()29|) . The final color trace in fundamental representation is trivial (and gives 1/2), and we derive 
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the relatively simple form 
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(31) 



in which all fields are in adjoint representation. There is only a single ordered exponential, linking 
the gauge index at the origin (c) with the field strength at the variable point A'. Note that the 
index c of $p is itself linked to the final state by the auxiliary ordered exponential that we have 
added in the I direction through gauge completion, as described above and in Ref. [T7j . 

To derive contributions of order v 2 , or equivalently q 2 /m 2 = 4q 2 /P 2 , we will study the nonlocal 
matrix element 
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(32) 



As above, (anti-) time ordering is implicit in the (complex conjugate) amplitudes. This is the 
complete 0(v 2 ) result for Z^ 1 ), Eq. (J2D, because, as we have seen above, at order O(v ), 
the quark and anitquark eikonal lines in both the amplitude or its complex conjugate cancel 
completely, and thus decouple for soft radiation. 

In the following, we will study the expansion of Eq. (}3*2*|) to NNLO. The explicit factor of A' in 
Eq. (}3*Tj) modifies the eikonal propagators. To see this, we can formally evaluate the A' integral in 
Eq. (}3~Tj) in terms of the Fourier transform of the field strength, F v ^ a {k) = J d 4 xF u ^ a (x) exp[— ik- 
x] . With this convention, momentum k flows into the field strength (and hence out of the eikonal 
lines). For a given order in the expansion of the adjoint ordered exponential in Eq. (J31j) . the 
lower limit of the A' integral is some value A m , the maximum value of A in the ordered fields 
p ■ A^(Xp) from &p A \ The relevant integral is then 
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(33) 



where we have integrated by parts. The second term gives a squared eikonal propagator. The 
first (boundary) term in brackets on the right-hand side gives the standard eikonal propagator 
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of Eq. (|26|) . times a factor of A m , producing a similar pattern in the next integral. The next 
integral will again give a squared propagator plus a boundary term, until the final A integral, for 
which the lower limit is zero and the boundary term vanishes. The result for a specific diagram 
is to replace the standard product of eikonal propagators by a sum of terms, in each of which 
one of the propagators is squared. Vertices for the operators p ■ A are unchanged. The relevant 
graphical notations for vertices are shown in Fig. [7| The three-point field strength vertex may 
be represented as 

tfiU(P> 9> fc ) = -9 5 - ((P • % M - (<? • k )P") > ( 34 ) 
and the four-point vertex as 

WZbdP, Q) = *9 2 Uc ( pV - W ) • (35) 

In both cases, c represents the color factor of the field strength tensor of Eq. 1)31)1. while a and/or 
b are the color indices of the gluon(s) that couple to the field strength. Because the adjoint 
eikonal lines end at the field strength in 1)31)1. corresponding to the color singlet pair in the final 
state, the three- and four-point vertices have only two and three color indices, respectively. 



P P Fqp P Fqp P Fgp 




Figure 7: a) Eikonal vertex; (b,c) vertices for the field strength; (d) line with eikonal and field 
strength interactions. 



As an example, corresponding to Fig. [7)1, we have the expression 
gP* U$JP, q, k) 1 1 ' 



-P ■ k' + ie (-P ■ k 



-P-k' + ie) 2 -P-k + ie 



(36) 



where we have chosen the sign of the infinitesimal imaginary part appropriate to the amplitude. 
To avoid clutter in and proliferation of figures, we will not introduce a graphical notation for 
squared propagators, but simply assume that the sum over terms is carried out in every diagram 
with a field strength operator at the largest time. 

The three-point vertex, Up in Eq. ()Hlj) is just the momentum representation of the Maxwell 
term of the field strength. It thus trivially decouples from scalar-polarized gluons, 



U F (p,q,k)k u = 0. 



(37) 



This result will lead to considerable simplification in our calculations; in particular, it eliminates, 
on a diagram-by-diagram basis, collinear poles associated with the octet eikonal line in the I 
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direction. This is because in any covariant gauge, collinear divergences are associated with 
gluons whose polarization is proportional to their momenta in the collinear limit |2*U] . 

In the next two sections, we apply these rules to study the coupling of soft gluons to the 
heavy quark pair. We begin at NLO, and then generalize to NNLO. 

5 Next-to-leading Order 

Figs. |Ht and b illustrate the origin of infrared divergences in the fragmentation function at 
next-to-leading order in a s to order v 2 . This infrared structure is the same as the lowest order 
contribution to X 2) Eq. (J32|) . As in Fig. El the sum over gluon connections to quark and antiquark 
on each side of the cut has been replaced by a single field-strength vertex. Because the parent 
gluon is off-shell by order m c , we may contract it to a point to study soft gluon corrections. For 
this purpose, it is then equivalent to study the matrix elements (J2"2j) to next-to-leading order, 
and that is how we shall describe our calculation below. We emphasize, however, that there is a 
trivial mapping from the matrix elements to the fragmentation functions. 

The vertical lines in Fig. |S] represent the quark-antiquark pair in the final state, and a pro- 
jection onto a color singlet (implemented by a color trace) is understood, along with a sum over 
all connections of the gluon to the quark and antiquark. The full set of diagrams is found by 
completing the cut, which can be done in only one way for|H^, where the gluon must be in the 
final state. For Fig. IHJd, on the other hand, there are two possibilities, one with a virtual gluon 
correction and one with a real gluon. In fact, of the two diagrams, onlylHk. can contribute to 
Eq. (j23). If we require a color singlet pair in the final state Fig. ISJd requires interference 
between octet and singlet in the hard scattering functions. We consider this diagram because it 
follows a pattern observed in the original arguments for NRQCD factorization, given in Ref. [7j, 
and because its square contributes to Z^ 1 ) at NNLO. 



Figure 8: Representative NLO contributions to g — ► cc fragmentation in eikonal approximation 
are found by all cuts of these diagrams. In these figures, the parent gluon is contracted to a 
point, represented by the dark circle, because it is off-shell by order m c . . 

Let us begin with Fig. |Ht, which has a topologically-factorized form, in which the soft gluon 
connects only to the heavy quarks, rather than to other finite-energy final-state lines. (In this 
case, the only such line is the eikonal line in direction /.) From the perturbative rules described 
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above, we immediately write down the following integral, which is readily evaluated in D = A — 2s 
dimensions, 



^ a \P,q) 



2,.2e 



16 g fi 



J^ZT^k 2 ) {q v {P-k)-{q.k)P v 

1 



x [q v {P -k)-{q- k)P u ] 



l(P-k) 2 } 2 



16 a s q 2 1 
3 vr P 2 -e 



(3? 



Here we have suppressed color factors, including the factor of 1/2 from the color trace mentioned 
above Eq. (|31|). The infrared pole in this result is familiar from NLO calculations of fragmentation 
|35| 13*6^ 13*7]. in which it is matched to the relevant NRQCD matrix element. 

We now turn to the cuts of Fig. |S1d, which contribute only to color interference terms. These 
diagrams, in which the gluon connects the quark-antiquark pair with the eikonal line, are not 
topologically factorized. Based on the arguments of |7j, we expect these to cancel, and they do. 
This was verified explicitly in Ref. [3B| f° r the case of color octet pairs in the final state. It will 
be instructive, however, to see how this happens in our velocity-expanded form to linear order in 
q with a color singlet final state, because the cancelation found here will be relevant to NNLO. 

We consider first the cut diagram with a virtual gluon loop in the amplitude. For our purposes, 
the overall normalization of the diagram is arbitrary, and we write 



,(86) 
J virtual 



with numerator factor 



d D k 
(2^p 



N(P,k,q, 



1 



(P-k + te) 2 k 2 + 



k + ie 



(39) 



N(P,k,q, 



2 [q ■ £(P ■ k) — P ■ £(q ■ k)] 
V2P r \V2k + q 3 + q ± -k ± 



(40) 



In this diagram, as in subsequent loop integrals, we will integrate first the minus loop momentum, 
by closing contours in the lower half-plane and picking up the relevant poles. Certain regularities 
and cancellations are conveniently represented in this manner, reducing the number of diagrams 
that must be computed explicitly. The result is shown in Fig. The double pole from the quark 
pair octet eikonal denominator, (P • k + ie) 2 , is always in the k~ lower half-plane, while the pole 
of the exchanged gluon is in the lower half-plane only when k + flows in the direction indicated 
in the figure. Closing in the lower half-plane, and neglecting the term odd in k± we find only a 
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terms proportional to 53, 

y (8b) _ gs g 2 

^virtual,! ~~ 1 ^ - 




00 dk + 



+ 



P (2tt) d - 1 J [J -00 (2k+ 2 + k\) 2 J 2k+ (k+ + k 2 ± /2k 



\2 



(41) 



These two terms correspond to the diagrams in Fig. |§k and b, where the straight line through 
the eikonal or gluon line indicates the k~ pole chosen. Both of these terms are logarithmically 
divergent by power-counting in the soft limit. As in the case above, however, they are collinear 
finite. In fact, the first term on the right-hand side is odd in k + and vanishes after symmetric 
integration. The pole at k + = 0, which would correspond to an on-shell intermediate state, has 
vanishing residue at order q. 

In the second term of the right-hand side of Eq. (|41j) . the exchanged gluon is on-shell with 
positive plus momentum flowing from the heavy quark pair to the eikonal lines. Its contribution 
to the cross section, as illustrated in Fig. Eb, is real, and it is straightforward to verify that it 
cancels the corresponding diagram for real gluon emission, illustrated by Fig. which has a 
standard flow of gluon momentum (k + > 0) in the cut diagram. We will use the notation shown 
in these figures to help organize our NNLO computations below. 

At the level of NLO, we have found that one-loop corrections indeed follow the expected 
pattern: they cancel except when topologically-factorized, and are thus consistent with matching 
to conventional NRQCD matrix elements. The presence of the octet Wilson line in our gauge- 
completed matrix elements does not change this pattern at NLO, as observed in Ref. [55] 

Before going on to the details of the NNLO calculations, we make a comment on gauge 
independence. As defined, the factorized fragmentation functions are gauge invariant, since the 
eikonal and the pair creation operators $^ ip^K, n x °f Eq. (|22|) are contracted to form a color 
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singlet vertex. Also, as we have seen, because of Eq. ()37|). there are immediate cancellations 
of many gauge terms like k^k u /(k 2 ) 2 for a gluon of momentum k, because of the field strength 
vertices that appear when we expand in the relative velocity of the quark pair. We can, of course, 
decouple the eikonal gauge line entirely, by choosing an i ■ A = gauge. The gauge invariance 
of the matrix elements assure that the result would be the same. For our purposes, however, 
Feynman gauge is most convenient. 

6 The Fragmentation Function at NNLO 

In this section we study in detail the infrared behavior of the gauge-completed gluon fragmenta- 
tion function at NNLO. Specifically, we will study non-topologically factorized diagrams at order 
a 2 in Eq. (|24|). We will find uncanceled infrared divergences for this set of diagrams, correspond- 
ing to 0(v 2 ) contributions to the fragmentation functions. We emphasize that the same infrared 
poles (proportional to oc 2 /(— e)) appear in cross sections, associated with soft gluon exchanges 
between the quark pair and a recoiling gluon. For this reason, the gauge-completion of matrix 
elements is necessary for factorization. At the same time, we will observe that the infrared pole 
is independent of the direction of the vector . This shows that the gauge-completed fragmenta- 
tion function is universal to NNLO. The same fragmentation function will match infrared poles 
for the quark pair recoiling against a gluon in any frame, or indeed (as we shall see), for any set 
of recoiling jets at this order of soft gluon exchange. We are not yet able to show, however, that 
this redefinition is universal at all orders in soft gluon exchange. 

6.1 The diagrams 

The soft-gluon diagrams that we will evaluate are shown in Fig. E3 Here, we compute only 
those contributions that correspond to the transition of a color octet pair to color singlet in both 
amplitude and complex conjugate. For NRQCD factorization to hold, all infrared divergences 
should either cancel or factorize into octet matrix elements. 

We will discuss the diagrams of Fig. El one at a time. In evaluating each diagram, k\ is 
defined as the momentum of the gluon that attaches to the quark pair at the left-most vertex, 
and it is always chosen to flow left to right in the diagrams. We label the momentum of the 
remaining gluon line attached to the quark pair as k%, choosing it to flow upward to the pair 
in each case. As observed in Sec. 5, the structure of the field strength vertex automatically 
eliminates collinear poles associated with soft gluons parallel to the i direction. 

For the purposes of this section, momenta associated with the quark eikonal lines, P, p = P/2 
and q, will all be scaled by the quark mass, m c . In the quarkonium rest frame, then, we have 
P = (2, 0), p = (1, 0) and q = (0, if/2). 

6.2 Summary of results 

In the remainder of this section, we have given the calculations that confirm our claims above 
in substantial detail. Since this discussion is of necessity rather detailed, it may be useful to 
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(I) 



(II) 



(III) 




+ c. c. 




+ c. c. 




+ c. c. 



(IV) 



(V) 



(VI) 



Figure 10: Diagrams I - VI discussed in the text. We sum over all cuts of these diagrams that 
can produce a color singlet quark pair. At the hard vertex, the pair is created in an octet state. 



summarize our results at the outset. Each of the diagrams in Fig. El contributes to the NNLO 
infrared factor in three separate quantities: first, the inclusive cross section for the production 
of a color-singlet charm pair of total momentum P to order q 2 in their relative momentum, 
through the fragmentation of an off-shell gluon; second, the fragmentation function for a gluon 
to the same color-singlet pair; and third, the gauge-completed production matrix element Eq. 
(122)1 for the production of the color-singlet pair from a local color-octet combination of quark 
and antiquark operators. The matching of the cross section with the fragmentation function was 
shown in Sec. 2, and the matching of the fragmentation function (fTK|) with the production matrix 
element (122)1 in Sec. 3. Notice, however, that these diagrams do not appear in the conventional 
production matrix element ([21)1 . Since all other diagrams are held in common between the two 
matrix elements (|22|) and ([21)1 . we can confidently conclude that the modification of the matrix 
element is necessary for matching at NNLO, at least in Feynman gauge. Indeed, because the two 
sets of diagrams are actually the same in a light-cone £ ■ A = gauge, this is a quick way to see 
that the matrix elements without the I eikonal lines are not gauge invariant. 

In the following subsections, we identify the diagrams by the numbers in Fig. EH I n view of the 
above, we recognize that each of the quantities: cross section, gluon fragmentation function, and 
production matrix element is proportional to the sum of these diagrams, multiplied by infrared- 
safe factors. There is no question that for individual final states, these diagrams are infrared 
sensitive. The question that we address in this calculation is whether, when all the final states of 
all the diagrams are combined, the infrared poles remain. As indicated above, the answer is yes. 

Discussing the diagrams one-by-one, the actual results (in covariant gauge) are rather simple 
to summarize. The infrared poles in dimensional regularization for diagrams I, II, IV, V and VI 
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all cancel. Only diagram III provides a noncancelling pole, given by 8 

2Re/// = -a 2 s 1 ^ , (42) 

with £f the relative velocity. Multiplied by the appropriate color factors, this pole will appear in 
the calculation of all of the three quantities just discussed. Eq. (J4*2~j) is the basic result of our 
calculation. Because it is nonzero, the gauge-completion of NRQCD matrix elements appears to 
be necessary to extend this formalism to production processes at NNLO. Conventional matrix 
elements simply will not match the infrared poles that are encountered in cross sections and 
fragmentation functions at this order. 

Having said this much, the reader who wishes to avoid, or delay, the details of the NNLO 
calulation may skip to the final subsection of this section, where we discuss how the applicability 
of this result to cross secctions involving the production of the pair with arbitrary numbers of 
hard jets, and to the conclusions, for a brief recapitulation. 



6.3 Ladder-like diagrams 

Diagrams I and II have a ladder and crossed-ladder structure. We discuss the calculation of 
infrared poles in II in some detail; diagram I has a very similar structure. We will show that the 
single IR pole of diagram II has an imaginary residule. 

The cuts of diagram II, that is the contributions from various final states, are shown in Fig. 
HTk . We begin with diagram IIA, in which a single gluon appears the final state. IIA is the 
complex conjugate of IIC, while IIB is real. Thus, at the order to which we work, we need 
consider only the real parts of each diagram. 

After dropping terms that are linear in ki±, i = 1, 2, the integral becomes 

[ (P ■ i) (g ■ h) - (q ■ t ) (P ■ h) ] [ (P ■ £) (g ■ k 2 ) - (g • I) (P • k 2 ) } 
(P • hY {P-k 2 - ie)\-t ■ k 2 - it) (I ■ (h - k 2 ) - it) 

_ r A dD r dD m 1 

x ^ -i ^ . (43) 

[k^ + k+) (k 2 + k+ - it) (-k+ + k+-ie) 

The first expression gives IIA in terms of P = (2, 0), normalized as in the perturbative expansion 
of Eq. (J3~2~j) . but suppressing color factors. In particular, as in the NLO case, we divide by (1/2) 2 
to compensate for the traces in quark representation. Here and below, A ~ m c is an ultraviolet 
cut-off for real soft gluon radiation. 

8 The factor of two on the left is a convention in our calculation below. 
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Again as in the one-loop examples, we do the minus integrals first, closing the k^ contour in 
the upper half-plane. This gives two terms, one from the quark-pair eikonal, another from the 
gluon propagator, 

II A = II A^ + II A^ , (44) 

represented in Fig. ITTb . Again as in the one-loop example, we readily verify that the gluon 
pole term, combined with the corresponding gluon pole of IIC, cancels the entire contribution 
of diagram IIB, where both gluons appear in the final state. For the remaining contribution we 
find after the minus integrations 



11 A ~ (2vr)2(^i) J d kl±d k2± J (2kt 2 + k 2 - 



2±> 

A A V. + 2 



4 k + 

x / dkt + ki ± m+ - k+ + ie ) ■ (45) 

In this expression it is clear that we may extend the ki± upper limit to infinity without changing 
the infrared behavior of the integrtal. Then, performing both transverse integrals, and changing 
variables to x = ^/2kf and y = k^/kf, we find 



IIA^ = 4(^) 2 (4vr) 2£ r 2 (l + £)g 3 2 



t,I ^L*«*(»-i+fc)' (46) 

This integral is infrared regularized for e < 0, that is in more than four dimensions. The pole is 
found from the identity 



1 1 £Y \ 

-6(x) + 



X l+Ne _ N£ 



+ 0(e) , (47) 



with a residue that is given by the e = limit of the remaining expression. The y integral has 
no pole, because for e < 0, the poles from y — > + and y — > _ cancel. The y integral at e = is 
then found to be 

Um "° L dy (v- 1 + (48) 

We conclude that although IIA^ k °\ and hence the complete diagram //, is infrared divergent, 
its divergence is imaginary, and does not contribute to the fragmentation function, which is real. 
Essentially identical considerations apply to the uncrossed ladder diagram, /. 



6.4 Diagrams with three gluons on the quark pair eikonal 

The diagrams with three gluons connected to the quark lines, are IV, V and VI of Fig. El We 
first consider diagram IV, which involves the commutator term of the field strength. Diagram V 
and VI both have an additional eikonal vertex at which a gluon couples to the quark pair in an 
octet color state. 



31 




Figure 11: (a) Cuts of diagram II; (b) k 2 poles of Fig. 2a. 



6.4.1 Diagram IV 





(IVB) 

Figure 12: Cuts of diagram IV. 



The relevant cuts of diagram IV are shown in Fig. As in the case of the ladder diagrams 
in the previous subsection, diagram IVB, which is real, cancels against the (k\ — k^) 2 = pole 
of IVA found by closing k^ in the upper half-plane. We thus need only evaluate the real part of 
IVA from the double pole at P ■ k 2 = 2k\ = 0. The imaginary part, of course, cancels against 
the complex conjugate diagram. 

In the same normalization as above, the numerator momentum factor for diagram IV is 
independent of k 2 , and is given by (after dropping terms linear in kx±), 



n IV = {P ■ k x q» - q ■ k^) (q^ - P^) F 
= A[-qlkt-{l/2)ql (fc+ + fcf)] , 



(49) 
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and the full contribution of diagram IV is given by the real part of 

IV A + IV B = IV A^ 

' a 



-32 



7T 



( 47r 2)l-2 E J Q ^+ 



dk2 



d D - 2 k 



d D - 2 k 



2_L 



-qlkt - (l/2)gi (kf + k\j2kt) 



(k+ + k\j2k\) 2 [-2(kt - k+)(kt + k\J{2kt)) - {k 21 - k 1A _f -ie] 2 ' 

(50) 

As in the previous diagrams, we change variables to y = k^/kf. In addition, we rescale the 
transverse integrals as /% = ki±/ \\/2k^). The Kj integrals are finite, but we find that an explicit 
infrared pole appears from the limit k± — > 0. After the k 2 integral, the diagram becomes 



IV = 4 



7T 



2e 3e-l 



X 



2 Z£ TT 



4 



dkt 

^.+ l+4e 



+ 



[l + «?] 2 2[1 + K ? 



«1 



1/(1 - k?) + + ie] 



l+e ' 



(51) 



The y integral is readily carried out (after changing variables to y' — y — (1/2) (1 — nf)), and 
we verify that the residue of the infrared single pole in e is imaginary. The real contribution of 
diagram IV to the fragmentation function is then infrared finite. 



6.4.2 Diagram V 

As for diagram IV, we will find the infrared pole of the real part of VA and VB, given in Fig. 
and once again the latter, with two gluons in the final state, cancels against the (k\ — k 2 ) 2 pole 
in the former, when the contour is closed in the upper half-plane. 





(VA) 

Figure 13: Cuts of diagram V. 



The Feynman rules for the field strength vertex lead to a sum of terms (in this case two) in 
which each of the denominators on the octet ordered exponential in Eq. (J31|) is squared, just as 
in Eq. (|36j). Thus, diagram VA has two terms. The momentum numerator factor, which depends 
only on k\ is the same for both. We choose to route the k\ momentum across the gluon to the 
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£-eikonal line on the bottom of the diagram, so that the right-most quark pair eikonal carries 
momentum k 2 to the right, and the exchanged gluon k 2 — ki up. The VA integral is then given 
by 



VA 



x 



i2 5 /y 

n VA (k 1 ,k 2 ) 
k-^ k 2 ^e 



A poo 
D 



d ki 



/OO 
d D k 2 5 + (kj: 



i 



(k 2 - ki) 2 - ie 



+ 



(fcf + kf) (kf + k^ - ie) (h{ + kf) (kf + k 2 - ie) 



(52) 



which exhibits the squaring of poles to the right of the field strength vertex in the diagram. The 
momentum numerator uva is 



n VA = (P • htf - q ■ hPV) (P • k igfl - q ■ k^) (P • £/2) 
= V2 ( -AqlktK - vl (H + K) 2 + 2(q ± ■ k 1± f ) . 



(53) 



After performing the k i integrals of VA and VB, and noting the cancellation of the exchange 
gluon pole, we are left with the contributions of the k\ pole 



V (k°) = 



dki 



+ roo 



dk+ I d°' 2 k 1± I d^k 



]D-2i 



( 47r 2)l-2e J q 2fc 4 

-2#1 - Ql (K + k\j2kf) 2 - 2(q ± ■ k 1± f 



(h + +klJ2kff 



x 



(2(kf -kf- ie)(-kf - kj ± /2ki) - (k 2± - k 1± ) 2 - ie) 
1 2 



{kf - ki - ie) (ki + k\j2kf) (2(kf - kf)(-kf - k{j2kf) - (k 2± - k 1± ) 2 - ie) 



(54) 

We rescale kf and both of the transverse momenta as y = kf /kf and Ki = ki±/V2kf, which 
again isolates an overall infrared divergence at the lower limit of the kf integration. The result 
can be expressed as 



*<* - W l£= f.$* f 



where the function J v is defined by 

/oo n 
dy / d 2 ~ 2e K 2 
-oo J 



4g|«| + qj (1 + Kp - Mg± ■ «ij 
(1 + k\Y 



(55) 



1 



(1 +K\){l-y- ie) {{y - l)(y + k\) + (k 2 - k,) 2 + ie) 
1 



[(y - l)(y + k\) + (k 2 - k^ 2 + ie] 



(56) 
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Performing the k 2 integral we find 



Jv(ki] 



-7T 



1-e 



r(i + £ ) 



dy 



[e (l + «f) (l _ j, - i e ) (j/2 + _ X ) 



(y 2 + - 1) - «? 



,1+e 



(57) 



The 1/e pole in the first term in square brackets comes from the term in which the squared 
denominator on the quark eikonal is outside the loop, and this pole is of ultraviolet origin. The 
corresponding k 2 virtual loop integral, as in the one-loop example of Fig. |HJd, Eq. ([39)1. is infrared 
finite, and hence may be absorbed into a coefficient function in the NRQCD expansion. At the 
same time, the remaining, k±, integral of this term is analogous to the k integral in Fig. |Sk, 
Eq. (|38|). and its infrared divergence is topologically factorized in the NRQCD expansion. Next, 
comparing the second term in brackets to Eq. (|5ip. we see that it is the same as the y integrand in 
that case, and when combined with the integral in Eq. (J53j) gives a purely imaginary infrared 
pole. In summary, the infrared sensitivity of diagram V is fully consistent with the NRQCD 
expansion. 

6.4.3 Diagram VI 

Diagram VI, with cuts shown in Fig. El is treated in a similar way to the previous two diagrams 
with three gluons connected to the eikonal quark pair line. In this case the k\ line is again 
connected to the left-most (field strength) vertex, while we route momentum k 2 — fci from the 
gluon eikonal (£) to the other field strength vertex. Once again the pole from the k 2 — k\ line of 
diagram VIA cancels the two-gluon final state, diagram VIB. 





Figure 14: Cuts of diagram VI. 
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The corresponding integral for VIA is 
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i2 5 / 2 g 4 
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ie 



(5? 



The first equality exhibits the squaring of poles to the right of the field strength vertex in the 
diagram. In the second, we note that the term in square brackets is a derivative with respect to 
kf, and that the expression is simplified by an integration by parts in that variable. 
The momentum numerator factor is 



n VIA = (P ■ krf - q ■ kxP") (Pj2) (P ■ (ki - k 2 )q v -q-(k x - k 2 )P v 
= 2V2q 2 3 (kf-k^)(kf-kf ) + ..., 



(59) 



where the terms linear in the ki± will not contribute, and are omitted in the second line. 

In the second form of Eq. (|58|). the k 2 integral has three simple poles. After the k 2 and 
integrals, using the cancelation of the exchange gluon pole, we have 
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(60) 

Carrying out our by-now standard rescalings and performing the k 2 transverse integration, we 
find 
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1-K 2 



i 



L (2(y - l) 2 - iey+z (2(y - l)(y + k\) - ie) 1 * 
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The two y integrals both give finite and purely imaginary contributions at e = 0, so that once 
again the sum of contributions to the fragmentation function from the cuts of diagram VI is 
infrared finite. We have now shown that of the six classes of diagrams generated from Fig. EH 
all but diagram III are consistent with standard NRQCD factorization. We now turn to this 
diagram, which is the most complex to compute. 

6.5 Three-gluon rescattering contribution 

Diagram III is distinguished by its three-gluon coupling. It connects a subdiagram analogous 
to Fig. |Hk, Eq. (jSHl)? which was infrared divergent but topologically factorized, with the eikonal 
line I. It describes a process in which the soft gluon that transforms the color octet pair to a 
color singlet pair rescatters on the adjoint eikonal to lowest order by exchanging a gluon. We 
recall that the gluon eikonal represents the influence of the remainder of the high-py process. 
We are thus testing the possible dynamical influence of this process on the soft hadronization 
itself. We shall find that it is a nontrivial influence, with a noncancelling infrared divergence. 
Nevertheless, the residue of the infrared poles will be rotationally invariant, and hence consistent 
with an NRQCD factorization in terms of our modified matrix elements. 
Before doing any integrals, diagram IIIA is of the form 



with a numerator factor njjj that we shall define below. As usual, we choose the rest frame of 
heavy quarkonium, P^ = (2, 0, 0, 0), and we will perform the fc7 integral by closing the contour 
in the upper half-plane. 

The basic pattern for diagram III in Fig. El is similar to those above: the two-gluon cut in 
Fig. EE IIIA, cancels the pole in from the exchanged gluon in IIIB that is attached to the 
octet eikonal line i. As for diagrams V and VI, we choose the momentum of this gluon as k\ — k2, 
flowing down. There are two additional poles in diagram IIIA when we close the k^ integral in 
the upper half-plane, as shown in Fig.EB After the cancellation with IIIB, only the contributions 
from poles (b) and (c) remain. 
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6.5.1 The numerator and the k\ pole 

The numerator factor nui is 

niii = (P ■ k x q» - q ■ kxP") (P • k 2 q v - q ■ k 2 P v ) 

x V [9nu{-h - k 2 ) a + gafi{2h - k 2 ) v + gva{2k 2 - h)n'. 
= 2 [-2(kt + k^)(q 2 k w k 20 - (q ■ hk 20 + q ■ k 2 k w )q + q ■ k x q ■ k 2 ) 
+4(7 • qk 20 - q ■ k 2 l°)(q ■ k 2 k 10 - q ■ kik 20 ) 

+4(7 ■ qk w - q ■ hl°)(q ■ hk 20 - q ■ k 2 k w )} 



This is a fairly complex expression, but is clearly symmetric in k\ and k 2 . 
When we take the contribution of the k\ pole, Fig. Hob , we find 

""*m - jj^ & > SM IM J >2 ^ J ' dD - % 

nm(k 1: k 2 ) k 2 =k 2 =0 

X (2(k+ - k+) (k 2 2± /2k+ - k\j2k+Y - (k 1± - k 2± ) 2 - ie) 



[kt + k 2 2± /2k+Y [kt + k\j2k\\ k i ~ k 2 
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This is an antisymmetric expression in k\ and k 2 , except for the imaginary contribution at 
kf = k 2 ■ As a result, once again all infrared poles in IIIA^ are imaginary and do not 
contribute to the fragmentation function. 



6.5.2 The k ( > = double pole 

We are left with the evaluation of the pole of diagram 116b. the double eikonal pole at k\ = 
as the only potential source of infrared singularities in the fragmentation function that are not 
topologically factorized in the usual sense. As we have anticipated, we will find an infrared pole 
in dimensional regularization. Since the calculation is a substantial one, we will give most of the 
details. To make it a bit more manageable, we first set q± =0 in the numerator (Jfillj) . and extend 
the result to nonzero transverse momentum in the appendix. 

At zero q±, the momentum numerator factor ()63|) simplifies to 

nm(«3,?i = 0) = 2q 2 3 [2(kt + k+)(kfk 2 + k+K) + A{k+ - k+){ktk 2 - k^k+)} . (65) 

We are now ready to pick up the pole in k 2 corresponding to the diagram of Fig. Mot , with the 
result 
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As above we work in D = 4 — 2e dimensions, and we rescale the transverse and k 2 momenta as 

h± k 2 ± 



«i 



V2kt 



2/ 



k 2 



+ ' 



(67) 
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which again isolates the infrared pole in the kf integral, 
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To do the k 2 integration we introduce a Feynman parametrization, 
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After the n 2 integration we get 
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We represent the above equation as 



IIIA^\q 3 ) = 2 2 +V- 1 ql (^)V £ 
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l+4e 



x J d 2 ^K X — [/«(«!) + /Wfa)], (71) 

where /W(fti) organizes a set of terms terms in which the pole at y = 1 has been cancelled by 
the numerator, and /^(fti) summarizes a set in which the factor remains. To effect this 
separation, we rewrite 1 — 3y = 3(1 — y) — 2 in the first term in brackets of Eq. ((70)) and to 
combine the second and third terms we use — x(l — y)y + (1 — x)y 2 = y(y — x). After these 
manipulations, we have 
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where we have completed the squares in the denominators. For the 1/(1 — y) terms we have 



J(2) ( 



«1, 



1-y-ie 



dx 



r(l + e) 



(y + x(k 2 - l)/2) 2 - x 2 {k\ + l) 2 /4 + ie ] 1+e 



T(2 + e) 



y(y-ar)(l - ft 2 ) 



(y + x(ft 2 - l)/2) 2 - x 2 (ft 2 + l) 2 /4 + it ] 2 + e 



(73) 



The y integral for /W, Eq. ([72)1. is straightforward. We change variables to j/' = y+x(n\ — 1)/2 
and note that in the numerator ?/(?/ — x) = y /2 + a; 2 (ftf — 1) — y'xnf, where the last term vanishes 
because it is odd in y' . In this way, we find 
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The x integration is now trivial, and using the expansion (— 1 +ie) £ 
isolate the imaginary pole in and a corresponding finite real part, 
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The real term in this expression, when substituted into Eq. (J71j) . gives a real, single pole, con- 
tribution to the fragmentation function from the kf integral. This is the generic mechanism we 
are after. 

The complete result at q± = 0, of course, requires I^(k{), which is a bit more complicated, 
because of the extra denominator 1/(1 — y — ie). We give the detailed calculation of 1^ in Sec. 
A.l of the appendix, where we show that 



/( 2 )( Kl ) = m - (1 -me) 



2k\ 



;i + k?)* 



+ .... (76) 



Substituting Eqs. ([75]) and (fTtj)) in Eq. (fTT]) we find that the remaining, K\, integration is conver- 
gent because K\ of order unity corresponds to ki± of order kf. Thus, the transverse momentum 
integration of the real gluon converges at a scale far below the fixed quark mass, and is effectively 
independent of the phase space cut-off. Performing the K\ integral, and adding the contribution 
from the complex conjugate of the diagram, we find for the leading, real - divergent term, 

2 Re IIIA^ P° le ) (q 3 ) = - a 2 s i qj (q ± = 0) . (77) 

The leading imaginary double pole, of course, cancels in the full fragmentation function. 

We evaluate the corresponding g^-dependent pole in the appendix. We note that all (73 x q± 
interference terms vanish because they are linear in the ki± integrations. Combining Eq. (J77J) for 
q± = with the result Eq. (jl(J5|) from the appendix for the q\ term, we obtain a rotationally 
invariant result 

2 Re IIIA^\q) = - a 2 s ±<f = - a 2 s 1 ^ . (78) 

This is the full result for diagram III and hence, as discussed in Sec. 6.2 above, for the entire 
NNLO infrared pole term in the cross section and fragmentation function, matched by the gauge- 
completed production matrix elements at the same order. 



6.6 Rotational invariance and universality 

The significance of rotational invariance is that the infrared pole is independent of the relative 
orientation of the pair's relative velocity v and the gluon eikonal direction I. The complete result 
shows first, that the gauge invariant redefinition of the NRQCD matrix element is necessary, 
but also shows that once this is done, the factorized form is consistent with universality of the 
factorization. 

As we have emphasized above, the same reasoning applies to cross sections in which the pair 
recoils against a gluon jet. In fact, the matching of cross sections with our matrix elements is 
even more general, as a result of the rotational invariance of Eq. ([78)1 . This follows from the 
nature of the gluon rescattering diagrams that give this result, in which two soft gluons attach 
to the pair, leaving only a single soft gluon to attach to the other jet. At the same time, the 
exchange of soft gluons at NNLO between the heavy quarkonium pair and each hard jet in the 
final state will give the same pole factor, given by Eq. (|78j). up to the effect of color. 
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Because the momentum factors are the same for the pole found by coupling the soft gluons 
to each final-state jet, we easily show that the complete color factor turns out to be independent 
of the number and directions of the jets, and of the color representations of their parent partons. 
We outline the proof as it applies to leptonic annihilation cross sections, where all jets are in the 
final state. The result follows from gauge invariance. 

Recalling the discussion of Sec. 2, we suppose that we are at a leading region of phase space 
where there is an arbitrary number of jets, of momenta pj, j — 1 . . . n. At any leading region, in 
the absence of soft gluon exchange, the cross section factorizes into a product of jet subdiagrams 
Jj(p 3)^,0,3, contracted in color indices aj and bj with a hard scattering function h ao , aj in the 
amplitude, and a corresponding function in the complex conjugate amplitude. In addition to the 
final-state jets, the short-distance function h is also contracted with the parent parton (gluon 
above) of the heavy quark pair, though color index a^. Concentrating just on the amplitude, the 
fragmentation function is thus proportional to the combination 



We suppress the function associated with the jet in which the pair appears. 

We now consider the effect of adding soft gluons at NNLO in this leading region, and we 
again discuss the case when the quark- ant iquark pair is an octet at short distances and a singlet 
in the final state. This requires that two gluons attach to the pair. Recalling the factorization 
property of jet-soft interactions derived in Sec. 2, the infrared behavior of each set of diagrams 
where soft gluons couple to jet j can be replaced by diagrams in which the soft gluons attach 
to an eikonal line in the direction of pj. Once again only soft-gluon diagrams like III in Fig. 
with a three-gluon coupling, can give rise to a real infrared pole in the cross section. In the set 
of such diagrams, the single exchanged gluon attaches to the n jets one at a time. 

We denote the color index of the exchanged gluon by e, and the flavor of the parent parton 
of jet j by fj. At fixed values of the pair relative velocity v = 2q, the effect of this insertion is to 
multiply M. of Eq. (|T§j) by the same pole term IIIA^ k °' pole \q), Eq. (J7HJ), that we have determined 
above for the fragmentation function at NNLO. This factor is independent of the jet to which 
the exchanged gluon attaches. 

The effect of the exchanged gluon's color, of course, differs from jet to jet, but is still quite 
simple after it has been factorized. The short- distance color tensor is multiplied by the matrix 
through which the soft gluon couples to the eikonal line in the pj direction, that is the color 
generator Tte>. In summary, the structure of the NNLO pole term in the cross section is 



We now observe that multiplication by the color generator for a given external line of the short- 
distance function h is equivalent to an infinitesimal color rotation of the corresponding external 
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line. The sum of color rotations on all its external lines vanishes by the gauge invariance of the 
theory. The sum of color rotations on all the final-state jets, therefore, is the negative of a color 
rotation on the parent gluon of the pair, 



n 



£5> 




a ,a ■ 



(81) 



1=1 a' 



The sum of the color factors associated with attaching a single soft gluon to all recoiling jets is 
therefore independent of the number and/or flavor of the final-state jets. The same argument 
can be applied to the color factors of the fragmentation function, with the same result. The 
gauge-completed matrix element is therefore universal up to NNLO, for arbitrary numbers of 
hard jets in the final state. 

7 Conclusions 

We have investigated the proposal of NRQCD factorization in production processes at large 
transverse momentumm (pt), and have demonstrated that factorization holds to NNLO in pro- 
duction from an octet pair, after a redefinition of the nonperturbative matrix elements in the 
effective theory. We have seen, in fact, that this matrix element is universal at NNLO for high-py 
quarkonium production with arbitrary final states. Many questions remain, however, and it is 
unclear to us whether the pattern we have found, uncanceled infrared divergences that can be 
absorbed into universal gauge-completed matrix elements, will survive at higher orders. On the 
other hand, the very nontrivial organization of the NNLO infrared divergences into a single power 
of q 2 is encouraging. 

So far our analysis has involved infrared structure associated only with electric dipole cou- 
plings, at momentum scales that are characteristically of order mv . Such a higher-order analysis 
will also require study of the lower momentum scale characteristic of binding energies mv 2 , which, 
as we have observed above do not enter into our NNLO octet-to-singlet calculations [33*] . Fi- 
nally, our study of fragmentation at large pt strongly suggests that the low-pr cross sections for 
quarkonium production deserves a fresh look [33] . 

In summary, the calculations and reasoning presented in this paper have, we believe, demon- 
strated that further investigation is crucial to provide a theoretical grounding for the analysis of 
the production of heavy quarkonia. 
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A Appendix 



The appendix provides more details on the evaluation of the infrared pole in diagram IIIA. The 
integrals presented here are all reasonably straightforward, and are complex only because of the 
rather large numbers of terms. Nevertheless, because they can be performed "by hand" , we feel 
an interested reader who wishes to reconstruct the calculation in detail may find the following 
relatively extensive presentation useful. 



A.l The integral 1^ of diagram IIIA with q± = 

We continue here with the detailed evaluation of the real a 2 e° contribution from the function 
J (2) (/ti) of Eq. (|Z3J). Compared to the case considered above, I^(ki), J^(kx) differs primarily by 
having an extra 1/(1 —y) denominator, which requires an additional Feynman parameterization. 

To simplify the y integral, we eliminate the explicit y 2 numerator factor by using the identity 
y 2 — yx = [y' 2 — x 2 (nf + l) 2 /4] — x(y — x)k 2 , where the term in brackets cancels a power in the 
denominator and where, as above, y' = y + x(k 2 — l)/2. This gives the slightly simpler form 
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where the second equality serves to define 1^ 2 '(k{) and j^{K\). Consider the first y integral of 
Eq. m, 
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Introducing an additional Feynman parameter, x', and expanding the square of the second de- 
nominator we get 
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The y integration is now easily performed, and gives 



r(i + e) v v 7o Jo *' 1/2 - £ 

X [x'x 2 ^ 2 + l) 2 /4 + (1 - x')\x{k\ - lj/2 + 1] + S^Hja^T ' 

To isolate the infrared pole of this expression, it is useful to change variables to u = x 2 and 
v' = (1 — x')/u. Also using (—1 + ie) = e l7r , we have 
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j(2) = _ 1 W^Tc; e -«r(3/2+ £ ) 



2T(l + £) Jo Jo (l-uv') 1/2 ' e 

X [(1 - ««')(«? + 1)74 + t/[X«? - 1)/2 + 1] + i(r^F^ ' 



The it and i>' integrals in (job)) are finite for e < 0, characteristic of an infrared pole. The 
1/e pole comes from w — > 0, and is isolated using (|47|). Its residue is purely imaginary. There 
is a corresponding real contribution to 1^ at £ — 0, however, found from the expansion of the 
exponential. The i>' integral is trivial at e = and u = 0, and we find 



Z< 2 > = m - (1 - me) 1 + 0(*e°) . (87) 

.1 + K-J 



This term will contribute at the level of a 2 /e in the fragmentation function after the integrals 
over K\ and kf. 

An identical procedure can be used to evaluate the second term, in Eq. ([FSjl. 

(k? + 1)/2 + V5«7[2(1-«i/)] 
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with the same denominator as in Eq. (J8b|) . The relevant singular behavior of this expression is 
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(89) 

Isolating the imaginary pole and the accompanying real finite part in the same way as for 
we get 
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Finally, substituting the results for and into Eq. (J82|) we find 
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which is the result quoted in Eq. (|7fi|). 
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A. 2 Transverse momenta for g in diagram III A 

When q 3 =0 the numerator for diagram III is 

nni(q 3 = 0,q ± ) = 2[ql{kt + kt){k+ + k^){kt + k 2 ) - 2{q L ■ k lA _) 2 {kt + k 2 ) 

-2(q± ■ k 2± ) 2 (kf + fcr) + 2( 9i . • • k 2 ±)(K + K)] ■ (92) 

Diagram IIIA^ with g 3 = is given by 
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Defining as above the scaled variables, Ki = , k 2 = and y = -^p, this integral becomes 
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which is analogous to Eq. (J58j) for the q 2 terms in the numerator. 

To perform k 2 integration we can again introduce a Feynman parameter x as in ([69)1 . The 
change of variables n 2 = k 2 — XK\ completes the square in the denominator, and the integral 
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results in the rather lengthy expression 
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where in the second relation, we define I^ 3 '(ki) to include the terms without the 1/(1 — y) 
denominator, and 1^ to include the remaining terms, all with this denominator. The infrared 
poles of I^ 3 '(ki) are identified in the same way as those of the corresponding q 2 integral, 
Eq. (JZ2J), while those of 7 (4) (ki) are found in the same way as for J (2) (/«i), Eq. (|ZSJ) - 
For I^ 3 ' the y integration is elementary, and we find 



J(3) ( 



«1, 



r(l/2 + e) 



r(3/2 + e; 



2 [-x 2 (l + /t 2 ) 2 + ie] 1 /2+ £ 

~(g± • ^i) 2 k 2 (l + «i) - - s 2 (k? - l)/2] 
[-a; 2 (l + /? 2 ) 2 + ie] 3 / 2 + £ 



(96) 



The overall infrared pole in this expression is easily identified as arising from the limit x —* 0. It 
comes entirely from the the middle term in the numerator of the second fraction, 



/ (3) («i) 



2m 



(q± ■ «i) 2 «i (l-i7re) 
(1 + kI 



2\3 



(97) 



Once again the pole is purely imaginary with, however, an associated finite real part. Finite 
corrections are all imaginary. 
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To evaluate I^(ki) we again introduce a Feynman parameter x', which enables us to do the 
K.2 integral just as for 1^ in Sec. A.l, giving 



/ (4) (^i) = 
T(2 + e) 



o jo 



dx I dx x' e 



dy 



(1 + Kl)ql - A(l - x)(q ± ■ K x y 



4 [x'[y 2 + xy(K 2 -l)-x 2 K 2 } + (l-x')(y-l) + ie] 2+e 



y(l — T(3 + e) -|— 



?/ (1 — x)xxT(3 + e) -j— 



[x'[y 2 + xy(n 2 - 


-1)- 


-x 2 k\] + (1- 


- x'){y - 


1) + ze] 2 + £ 




(q±- 


Kl) 2 [x(l + K 2 ) 


- «?] 




[x'[y 2 + xy(n 2 


-1) 


- X 2 /? 2 ] + (1 - 


- 


1) + ie] 3+£ 














[x'[y 2 + xy(nl 


-1) 


-x 2 k\] + (1 




- 1) + ze] 3 + £ 



(9? 



These integrals are precisely of the form of those in I^ 2 \ but to limit the rather large number of 
terms, we introduce the y-independent quantities 

P = \ [*(«?- 1) + (l-x')/x'} 

M 2 = x 2 k\ + {l-x')/x. (99) 
In this notation, the denominators of Eq. ()98|) are 

x'[y 2 + xy{n\ - 1) - x 2 k\} + (1 - x'){y - 1) = x(y 2 + 2Py - M 2 ) . (100) 
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After the y integral, I^(ki) can be written as a sum of five terms, 



I (4) (^i) = • 
T(3/2 + e 



irte 



x 



dx I dx' 
o Jo 



1 + K\)ql 



x 



[pa + M 2 } 3 ^ 



r (3/2 + e) (l + «j)gj 

4 1 J x' 2 [P 2 + M 2 ] 3 /2+ £ 



-r(3/2 + e) (1 - x) - 



-r(5/2 + e)x(l -x)P 

r(3/2 + £ 



x /2[p2 + M 2]3/2+ £ 

(g± ■ K!) 2 (x - (l -x)k 2 ) 



x(l — x) 



x /2[p2 + M 2]5/2+ £ 

(g ± -K 1 ) 2 (M 2 -2P 2 + eM 2 ] 



= — \/-Kie 



r{3/2 + e)[i®(K 1 ) + + iW( Kl ) + z w (ki) + z^(ki)], 



x /2[p 2 + M 2]5/2H 
,'( 2 ) 



(3), 



,•(5), 



(101) 



where the final line is the notation we will use for the five terms, taken in order, with v& the 
first, and the last. 

The infrared pole of each of the can be found by the straightforward, if slightly tedious, 
application of the following steps: 1) re-express P and M 2 in terms of x, x' and k\ using ()99j) . 2) 
change variables as above to u = x 2 and v' = (l—x')/u, 3) identify the residue of the singular 
u~ l ~ s behavior, where it is present. In fact, of the five terms, only and are singular 

at e = 0. Their poles are determined from 



du 

o u JO 



dv' 



[v> + (1 + ^) 2 /4] 3 / 2 ' 



(102) 



^ (3) (^; 



du 



u 



l+e 



dv' 



[v> + (1 + K 2 ) 2 /4] 3 / 2 ' 



(103) 



* {4) k: 



- 2 (4 - 1) 



du 



u 



l+e 



dv' 



r 2\2 



/4]5/2 



(104) 



The remaining two terms behave as u~ l l 2 for u — > and are hence of order e° and imaginary. 

The contributions of the poles in Eqs. (jl02|) - (|104p are also purely imaginary because of the 
overall factor of —in in (jlOlj) . and will enter the fragmentation function as an imaginary double 
pole, which therefore cancel. Correspondingly, all real terms of order e° from the u, v' and k± 
integrals contribute only at the level of an imaginary single pole. A real single pole in the final 
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result can only result from a relative factor —is, which is found as above from the expansion of 

g— ins 

The final result for IIIA^ k °\q±), defined by Eq. (}9~3*j) and (}9~2*j) . is therefore found from: 4) 
isolating the finite real part from the expansion of the overall factor of e~ l7TS in (jlOij) . 5) performing 
the remaining K\ integration at e — 0, and finally 6) replacing the final kf integral by 4e), 
according to Eq. ([47)1 . In this way, we obtain 

2 Re IIIA^° le \q ± ) = - a] 1 £ («& = 0) , (105) 

matching Eq. (|77|) for the g| term. There are no terms linear in and and the complete 
result is thus rotationally invariant. 



References 

[1] A. De Rujula, in 50 years of Yang-Mills theory*, ed. G. 't Hooft (World Scientific, Singapore, 
2004) p. 401, |hep-ph/0404215l 

[2] E. Braaten and T.C. Yuan, Phys. Rev. Lett. 71, 1673 (1993); Phys. Rev. D50, 3176 (1994); 
E. Braaten, K. Cheung and T. C. Yuan, Phys. Rev. D48, 4230 (1993); Y.Q. Chen, Phys. Rev. 
D48, 5181 (1993); T.C. Yuan, Phys. Rev. D50, 5664 (1994); E. Braaten, M.A. Doncheski, 
S. Fleming and M. Mangano, Phys. Lett. B333, 548 (1994); D.P. Roy and K. Sridhar, Phys. 
Lett. B339, 141 (1994); M. Cacciari and M. Greco, Phys. Rev. Lett. 73, 1586 (1994). 

[3] A.H. Mueller, Phys. Rev. D18, 3705 (1978). 

[4] J.C. Collins and G. Sterman, Nucl. Phys. B185, 172 (1981). 

[5] J.C. Collins, D.E. Soper and G. Sterman, in Perturbative quantum chromodynamics, ed. 
A.H. Mueller (World Scientific, Singapore, 1989), p. 1; | |hep-ph/0409313| 

[6] J.C. Collins, D.E. Soper and G. Sterman, Nucl. Phys. B261, 104 (1985); B308, 833 (1988); 
G.T. Bodwin, Phys. Rev. D31, 2616 (1985) Erratum ibid. D34, 3932 (1986). 

[7] G.T. Bodwin, E. Braaten and G. P. Lepage, Phys. Rev. D51, 1125 (1995), Erratum ibid. 
D55, 5853 (1997) |hep-ph/9407339| 

[8] E. Braaten, S. Fleming and T.C. Yuan, Ann. Rev. Nucl. Part. Sci. 46, 197 (1996) 
|hep-ph/9602374| 

[9] N. Brambilla et al. (Quarkonium Working Group) hep-ph/0412158, and references therein. 

[10] E. Braaten and S. Fleming, Phys. Rev. Lett. 74, 3327 (1995); E. Braaten, S. Fleming and 
T. C. Yuan, Ann. Rev. Nucl. Part. Sci. 46, 197 (1996); P. L. Cho and A. K. Leibovich, 
Phys. Rev. D53, 6203 (1996); Phys. Rev. D53, 150 (1996); E. Braaten, S. Fleming and A. 
K. Leibovich, Phys. Rev. D63, 094006 (2001). 



51 



[11] E. Braaten, B. A. Kniehl and J. Lee, Phys. Rev. D62, 094005 (2000); E. Braaten and J. Lee, 
Phys. Rev. D63, 071501 (2001); M. Beneke and M. Kramer, Phys. Rev. D55, 5269 (1997); 

A. K. Leibovich, Phys. Rev. D56, 4412 (1997). 

[12] M. Cacciari and M. Kramer, Phys. Rev. Lett. 76, 4128 (1996); M. Beneke, M. Kramer and 
M. Vanttinen, Phys. Rev. D57, 4258 (1998); J. Amundson, S. Fleming and I. Maksymyk, 
Phys. Rev. D56, 5844 (1997); R. M. Goodbole, D. P. Roy and K. Sridhar, Phys. Lett. B373, 
328 (1996); B. A. Kniehl and G. Kramer, Phys. Rev. D56, 5820 (1997). 

[13] C. G. Boyd, A. K. Leibovich and I. Z. Rothstein, Phys. Rev. D59, 054016 (1999); M. Klasen, 

B. A. Kniehl, L. N. Mihaila and M. Steinhauser, Phys. Rev. Lett. 89, 032001 (2002). 

[14] G. C. Nayak, M.X. Liu and F. Cooper, Phys. Rev. D68, 034003 (2003); F. Cooper, M.X. 
Liu and G. C. Nayak, Phys. Rev. Lett. 93, 171801 (2004); M. Klasen, B. A. Kniehl, L. N. 
Mihaila and M. Steinhauser, Phys. Rev. D68, 034017 (2003), 

[15] M. Beneke and I. Z. Rothstein, Phys. Rev. D54, 2005 (1996) [Erratum-ibid. D54, 7082] 
(1996)]; W. K. Tang and M. Vanttinen, Phys. Rev. D54, 4349 (1996); S. Gupta and K. 
Sridhar, Phys. Rev. D54, 5545 (1996). 

[16] M. Kramer, Prog. Part. Nucl. Phys. 47, 141 (2001) |hep-ph/0106"T20l 

G.T. Bodwin, at 2nd International Conference on Flavor Physics (ICFP 2003), Seoul, Korea, 



6-11 Oct 2003, hep-ph/0312173 



[17] G.C. Nayak, J.-W. Qiu and G. Sterman, Phys. Lett. B 613, 45 (2005) |hep^ph/0501235| 
[18] A. Czarnecki and K. Melnikov, Phys. Rev. Lett. 80, 2531 (1998) |hep-ph/9712"222l 
[19] M. Beneke, A. Signer and V.A. Smirnov, Phys. Rev. Lett. 80, 2535 (1998) |hep^pli79 712302 
[20] G. Sterman, Phys. Rev. D17, 2773 (1978), ibid. 2789 (1978). 

[21] G. Sterman, An introduction to quantum field theory (Cambridge Univ. Press, Cambridge, 
1993), chap. 13. 

[22] S. Coleman and R.E. Norton, Nuovo Cim. 38, 438 (1965). 

[23] R.J. Eden, P.V. Landshoff, D.I. Olive and J.C. Polkinghorne, The analytic S-matrix (Cam- 
bridge Univ. Press, Cambridge, 1966), chap. 2. 

[24] J.M.F. Labastida and G. Sterman, Nucl. Phys. B254, 425 (1985); J.W. Qiu and G. Sterman, 
Nucl. Phys. B353, 137, (1991). 

[25] J.C. Collins and D.E. Soper, Nucl. Phys. B193, 381 (1981). 

[26] C.W. Bauer, D. Pirjol and I.W. Stewart, Phys. Rev. D65, 054022 (2002) |hep-ph/0109"045l 



52 



[27] C.W. Bauer, D. Pirjol and I.W. Stewart, Phys. Rev. D66, 054005 (2002) |hep^h70 205289 

[28] R. Basu, A.J. Ramalho and G. Sterman, Nucl. Phys. B244, 221 (1984). 

[29] G.T. Bodwin, S.J. Brodsky and G.P Lepage, Phys. Rev. Lett. 47, 1799 (1981). 

[30] S.B. Libby and G. Sterman, Phys. Rev. D18, 3252 (1978). 

[31] J.C.Collins, L. Frankfurt and M. Strikman, Phys. Lett. B307, 161 (1993) |hep^h/ 9212212} 
E.L.Berger, J.C. Collins, D.E. Soper and G. Sterman, Nucl. Phys. B286, 704 (1987). 

[32] C.W. Bauer, D. Pirjol and I.W. Stewart, Phys. Rev. D68, 034021 (2003) |hep-ph/0303T56"l 

[33] J.C. Collins and D.E. Soper, Nucl. Phys. B194, 445 (1982). 

[34] J.C. Collins and W.-K. Tung, in proceedings of the 23rd International Conference on High- 
Energy Physics, Berkeley, CA, 16-23 Jul y 1986, N ucl. Phys. B278, 934 (1986); J.C. Collins 
Phys. Rev. D58, 0940021 (1998) |hep^ph/9806259[ 

[35] J.P Ma. Phys . Lett. B332, 398 (1994) hep-ph/940124; Nucl. Phys. B447, 405 (1995) 
|hep-ph/9503346"l 

[36] E. Braaten and Y.-Q. Chen, Phys. Rev. D55, 7152 (1997) |hep^h/9701242[ 

[37] A. Pet relli, M. Cac ciari, M. Greco, F. Maltoni and M.L. Mangano, Nucl. Phys. B514, 245 
(1998) |hep^p"h/9707223| 

[38] J. P. Ma and Z. G. Si, |arXiv:hep-ph/0506078| 

[39] E. Braaten and J. Lee, Nucl. Phys. B586, 427 (2000) |hep-ph/0004"228l 
[40] J. Lee, Phys. Rev. D71, 094007 (2005) |hep^ph/0504285[ 

[41] M. E. Luke, A. V. M anohar and I. Z. Rothstein, Phys. Rev. D 61, 074025 (2000) 
|arXiv:hep-ph/99102"09| . 

[42] E. Braaten, lectures at the Third International Workshop on Particle Physics Phenomenol- 
ogy, Taipei, November 1996, |hep-p h/9702225, 

[43] N. Brambill a, A. Pineda (CERN), J. Soto and A. Vairo, Nucl. Phys. B566, 275 (2000) 
|hep-ph/9907240| 

[44] V.A. Khoze, A.D. Martin, M.G. Ryskin and W.J. Stirling, Eur. Phys. J. C39, 163 (2005) 
hep-ph/041002. 



53 



